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Abstract 


A  computer  model  of  a  pressurized  water  reactor  (PWR)  was  developed  for  use  as  a 
teaching  tool  in  graduate  level  nuclear  reactor  courses.  The  development,  based  on  the 
diffusion  equation,  includes  the  methodology  for  solving  the  steady  state  spatial 
dependence  of  the  neutron  power  output  in  a  homogeneous  right  circular  cylinder 
unreflected  PWR  system.  This  includes  a  two  dimensional  one  energy  group  model,  a 
three  dimensional  one  energy  group  model,  and  a  three  dimensional  two  energy  group 
model. 

To  solve  the  homogeneous  diffusion  equation,  a  method  was  developed  to  search  for 
criticality  of  the  reactor  based  on  the  geometry  and  reactor  core  material  composition. 
For  the  one  energy  group  models,  a  perturbation  technique  was  developed  to  assist  the 
program  user  in  modifying  the  macroscopic  absorption  coefficient  to  drive  the  reactor  to 
criticality.  For  the  three  dimensional  models,  a  blocked  tridiagonal  solver  was 
developed  to  solve  the  numerical  linear  system  of  equations  approximating  the  diffusion 
equation. 

The  model  was  coded  using  Visual  BASIC  5.0^'^.  This  provides  a  platform  that  is 
exportable  to  personal  computers  and  allows  direct  linkage  to  Windows  based  programs. 
The  code  automatically  charts  and  displays  the  power  distribution  profile  using  ExceF'^ 
and  displays  the  calculated  multiplication  factor  determining  criticality. 
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MODELING  PRESSURIZED  WATER  REACTOR  KINETICS 


I.  Introduction 


Background 

Basic  nuclear  reactor  courses  at  the  graduate  and  undergraduate  level  focus  on 
teaching  students  how  to  calculate  radial  and  axial  flux  and  power  for  steady  state  (non¬ 
time  dependent)  reactors.  Many  nuclear  reactor  textbooks  cover  the  fundamentals  of 
nuclear  physics  and  apply  the  diffusion  equation  to  approximate  the  behavior  of  neutrons 
within  the  reactor  core.  Typically,  the  reactor  books  review  one-dimensional,  one  speed, 
homogeneous  models  for  various  geometric  shapes  in  great  detail.  Some  even  outline 
numerical  approaches  for  solving  the  approximate  solution. 

While  students  often  gain  a  basic  understanding  of  the  general  physics,  they  typically 
lack  a  qualitative  and  intuitive  understanding  of  the  reactor  core  nucleonic  behavior  based 
on  core  geometry  and  composition.  A  computer  model  can  be  used  to  provide  the 
students  with  a  tool  that  can  visually  explain  how  flux  and  power  are  impacted  by 
changing  the  core  geometry  and  composition. 

Problem  Statement 

The  problem  statement  for  this  thesis  was  to  develop  a  working  computer  model  of  a 
pressurized  water  reactor  (PWR)  for  use  as  a  teaching  tool  in  graduate  level  nuclear 
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reactor  courses.  The  development  includes  the  equations  and  methodology  for  solving 
the  steady  state  spatial  dependence  of  the  neutron  flux  and  power  output  in  a 
homogeneous  right  circular  cylinder  unreflected  PWR  system.  This  includes  a  two 
dimensional  one  energy  group  model,  a  three  dimensional  one  energy  group  model,  and  a 
three  dimensional  two  energy  group  model. 

Approach 

The  fundamental  approach  was  to  model  the  reactor  using  the  diffusion  equation.  For 
a  steady  state  system,  the  diffusion  equation  reduces  to  the  Helmholtz  equation. 

VV(r)  +  B,V(r)  =  0 

where 

BI  =  Geometric  buckling  =  ^ ^  (1/cm^). 

*  D 

Since  this  is  a  homogeneous  equation,  one  must  determine  the  eigenvalues  to  achieve  a 
non-trivial  solution.  For  cylindrical  geometries,  the  eigenfunction  corresponding  to  the 
smallest  eigenvalue  is  non-negative  everywhere  within  the  reactor.  This  is  physically  the 
only  value  of  importance  because  the  flux  cannot  be  negative  in  a  reactor. 

To  achieve  a  physical  solution,  we  rewrite  the  Helmholtz  equation  and  insert  an 

eigenvalue  . 

-V.Z)V<Zi(r)  +  <Pir)  =  ^ V  S, 

Where 

V  =  average  number  of  neutrons  per  fission 

Sa  =  the  probability  of  absorption  per  unit  path  length  (1/cm). 
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For  a  particular  value  of  k,  this  equation  will  have  a  unique  solution.  As  will  be  shown 
later,  if  k  equals  one  the  reactor  is  critical.  If  k  does  not  equal  one,  the  core  geometry 
and/or  material  composition  must  be  changed.  Searching  for  the  flux  when  k  equals  one 
is  called  the  criticality  search.  This  criticality  search,  using  the  diffusion  equation,  is  the 
basis  for  the  development  of  the  code. 

The  first  step  in  developing  the  code  was  to  solve  the  two  dimensional,  one  energy 
group  diffusion  equation  using  the  finite  central  difference  method.  The  finite  central 
difference  method  provides  a  satisfactory  order  of  accuracy  and  is  generally  used  as  the 
initial  method  for  modeling  or  designing  reactors.  The  finite  difference  method  results  in 
solving  a  tridiagonal  matrix  system  using  a  power  iterative  technique  to  solve  for  the  flux 
at  criticality.  The  program  in  this  thesis  uses  the  Crout  factorization  method  to  solve  the 
tridiagonal  system  of  equations.  A  perturbation  technique  is  used  to  perturb  an  initial 
guess  of  the  macroscopic  cross  section  to  drive  the  modeled  reactor  to  a  critical  level. 

This  perturbation  will  assist  the  user  in  selecting  the  macroscopic  cross  section  that  will 
result  in  a  critical  condition.  The  one  group  model  assumes  that  the  energy  of  the 
neutrons  is  equal  at  every  spatial  point  within  the  reactor.  The  model,  based  on  a 
homogeneous  un-reflected  reactor,  which  is  not  time  dependent,  yields  a  two  dimensional 
solution  of  power  versus  radial  position.  The  initial  two  dimensional  model  was 
developed  using  FORTRAbF^  and  then  later  converted  to  Visual  BASIC  5.0 

The  two  dimensional  one  group  model  was  expanded  to  a  third  dimension  by  adding 
a  solution  in  the  axial  direction.  The  output  of  this  model  provides  both  radial  and  axial 
power  plots  in  three  dimensions.  I  used  the  finite  central  difference  method  to 
approximate  the  diffusion  equation.  This  changed  the  system  from  a  pure  tridiagonal  to  a 
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blocked  tridiagonal  system  because  of  the  additional  sub  and  super  diagonals.  This 
required  the  development  of  a  blocked  tridiagonal  solver  to  solve  the  system  of  equations 
and  provide  the  flux  at  each  interior  mesh  point.  With  all  flux  values  known  at  the 
interior  points,  both  the  axial  and  radial  power  distributions  can  then  be  plotted.  The 
three  dimensional  model  was  written  in  Visual  BASIC  5.0. 

The  final  step  was  to  develop  a  two  energy  group  three  dimensional  model.  The 
model  assumes  only  down  scatter  of  neutrons  that  are  directly  coupled,  meaning  neutrons 
only  scatter  to  the  next  lowest  energy  level.  This  model  uses  the  finite  central  difference 
method  and  the  blocked  tridiagonal  solver. 

Visual  BASIC  5.0  was  chosen  because  it  offers  many  advantages  over  scientific 
languages  such  as  FORTRAN.  It  allows  the  programmer  to  build  an  executable  file  that 
links  automatically  into  simple  plotting  tools  such  as  Excel.  With  Visual  BASIC  5.0,  you 
can  command  and  control  Excel  as  well  as  other  Microsoft  Windows  software.  For 
example,  the  Visual  BASIC  5.0  reactor  code  populates  an  Excel  spreadsheet  with  the 
solution  data  and  then  builds  the  charts  all  from  within  Visual  BASIC  5.0.  The  charts  are 
linked  and  updated  to  appear  as  an  object  on  the  Visual  BASIC  5.0  form.  This  capability 
provides  the  user  with  automated  graphs  of  the  power  based  upon  the  input  parameters  as 
well  as  access  to  the  output  data  on  spreadsheets.  This  advantage  precludes  the  user  from 
having  to  manually  create  the  charts  or  plots  in  another  computer  language.  These 
features  outweigh  the  advantages  of  FORTRAN  such  as  computational  speed  and  built-in 
intrinsic  functions.  Additionally,  one  can  export  the  program  packaged  with  the  runtime 
dynamic  link  language,  thus  not  requiring  Visual  BASIC. 
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The  advancements  in  computer  technology  have  made  using  Visual  BASIC  5.0  an 
alternative  to  scientific  programs.  Less  than  five  years  ago  personal  computers  were  too 
slow  to  solve  three  dimensional  diffusion  problems  using  Visual  BASIC  5.0.  The  low 
cost  and  improvements  of  memory  and  processors  allow  personal  computers  to  be 
capable  of  solving  complex  numerical  problems  using  Visual  BASIC  5.0  in  a  fraction  of 
previous  times. 
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II.  Theoretical  Development 


Diffusion  Theory 

“Reactor  kinetics  is  the  area  of  reactor  physics  concerned  with  predicting  what 
happens  to  the  neutron  flux  density  when  the  balance  condition  associated  with  the 
critical  state  is  disturbed  (Henry,  1986:296).”  The  generation  of  heat  in  a  reactor  system 
is  proportional  to  the  fission  rate,  which  is  a  function  of  the  neutron  flux.  The  neutrons  in 
a  thermal  reactor  range  in  energies  from  0.001  eV  to  about  10  MeV.  To  simplify  the 
design  process  of  reactors,  neutrons  are  divided  into  energy  groups.  The  one  group 
model  deals  with  the  thermal  neutrons  only;  however,  it  also  accounts  for  those  produced 
from  both  prompt  and  delayed  neutrons.  The  two-group  model  deals  separately  with  both 
thermal  and  fast  neutrons. 

It  is  common  practice  to  approximate  the  exact  neutron  transport  equation  using 
diffusion  theory.  The  neutron  transport  equation  accounts  for  the  angular  dependent 
neutron  density  within  a  volume.  The  diffusion  equation  is  the  result  of  removing  the 
angular  dependence  from  the  transport  equation. 

The  diffusion  equation  is  based  on  Pick’s  Law  and  the  equation  of  continuity.  Pick’s 
law  is  shown  in  equation  (1). 


J.= 


(1) 


where 

=  the  net  number  of  neutrons  passing  a  unit  area 
perpendicular  to  the  x-direction  in  a  unit  of  time 
D  =  the  diffusion  coefficient  (cm) 

(j)  =  the  flux  (neutrons/cm^)(cm/sec) 
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Pick’s  law  was  originally  used  to  predict  the  flow  of  chemicals  from  one  region  of 
higher  concentration  to  another  region  of  lower  concentration  solute.  The  flow  was 
found  to  be  equal  to  the  negative  gradient  of  the  solute  concentration.  Although  neutrons 
do  not  actually  flow,  their  behavior  can  be  modeled  using  this  concept  (Lamarsh, 

1983: 192).  Early  reactors  were  designed  using  this  technique.  Today,  more 
sophisticated  and  computationally  demanding  methods  are  available  to  design  reactor 


cores. 


To  develop  the  diffusion  equation  one  begins  by  using  the  equation  of  continuity. 
The  equation  of  continuity  states  that: 


The  rate  of  change  in 
number  of  neutrons  per 
volume  (V) 


production  rate 
of  neutrons  in  V 


absorption  rate 
of  neutrons  in  V 


leakage  rate  of 
neutrons  in  V 


.(2) 


By  substituting  Pick’s  law  into  the  equation  of  continuity,  the  general  diffusion  equation 
becomes: 


+  (!>  =—,  (3) 

^  dt 

where  equation  (3)  is  the  non-steady  state  diffusion  equation  and 

D=  the  diffusion  coefficient  (cm) 

=  the  Laplacian  (divergence  of  the  gradient) 

^  =  the  neutron  flux  (neutron  cm/cm^  sec) 

Ea  =  macroscopic  absorption  cross  section  (1/cm) 

Ef  =  macroscopic  fission  cross  section  (1/cm) 

V  =  neutrons/fission. 

Removing  the  time  dependence  results  in  the  Helmholtz  equation. 
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-DVV  +  S>  =  vX/^  (5) 

This  is  the  fundamental  equation  to  be  solved  for  the  solutions  to  the  problem  statement. 

The  development  of  the  boundary  conditions  is  key  to  the  solution  of  the  diffusion 
equation  for  finite  cylindrically  shaped  reactor  cores.  In  order  to  develop  a  physical 
meaning,  the  total  flux  must  be  positive  and  real  in  all  areas  within  the  core.  The 
diffusion  equation  and  Pick’s  law  are  not  valid  at  physical  boundaries  since  they 
approximate  the  value  several  mean  free  paths  inside  the  boundary.  To  account  for  the 
physical  boundaries,  the  diffusion  method  models  the  measured  flux  by  assuming  the  flux 
is  zero  at  an  extrapolated  distance  beyond  the  outer  physical  boundary  layer  of  the  reactor 
core.  The  exact  flux  does  not  reduce  to  zero  beyond  the  boundary;  however,  the 
diffusion  theory  assumption  allows  for  reasonable  flux  calculations  within  a  few  mean 
free  paths  of  the  boundary  (Duderstadt  and  Hamilton,  1976:144).  See  Figure  1  for  a 
graphical  comparison  of  the  measured  flux  and  the  diffusion  theory. 

The  extrapolated  distance  for  plane  geometries  is  calculated  by  using  equation  (6). 

rf- 0.714  (cm) 


where  the  transport  mean  free  path  is 


4=3D=  — 


8 


Figure  1  Extrapolated  distance  at  outer  boundary 

However,  for  relatively  large  reactors  the  extrapolated  distance  can  be  neglected 
without  significantly  impacting  the  order  of  accuracy  because  the  extrapolated  distance  is 
on  the  order  of  centimeters  or  less  as  compared  to  the  radius  of  approximately  one  to  two 
meters  on  average.  In  this  model,  the  assumption  is  made  that  the  flux  is  zero  at  the  top, 
bottom,  and  sides  of  the  reactor  core  and  the  derivative  of  the  axial  and  radial  flux  at  the 
centerline  of  the  reactor  equals  zero.  This  is  accomplished  by  setting  the  flux  at  the 
centerline  equal  to  the  flux  at  the  first  interior  mesh  point  away  from  the  centerline. 

The  multi-group  diffusion  equation  discretizes  the  range  of  neutron  energies  into 
energy  groups  as  shown  in  Figure  2.  Notice  that  the  grouping  begins  with  the  highest 
energy  group  number  and  works  toward  the  lowest  energy  group  number.  The  highest 
energy  group  number  corresponds  to  the  lowest  energy  level  of  the  neutrons. 
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Figure  2  Energy  Groups 


Equation  (7)  can  be  described  by  the  energy  dependent  version  of  the  diffusion 
equation.  The  equation  is  based  on  integrating  the  neutron  energy  (averaging)  over  the 
energy  group  of  concern,  Eg<E<  Eg.i. 


The  rate  of  change 
of  neutrons  in 
Group  g 


Change  due  to 
leakage 


absorption  in 
group  g 


source 
neutrons 
in  group  g 


neutrons 
scattering 
out  of  group  g 


+ 


neutrons 
scattering  into 
group  g 


(7) 


For  the  two  energy  group  model,  the  energy  groups  are  shown  in  Figure  3. 
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Thermal 


Fast 


g=2 


g=l 


E2=0eV  E,=leV  Eo=10MeV 


Figure  3  Energy  Spectra  for  Two  Energy  Group  Model 


The  development  of  the  two  energy  group  system  is  based  on  the  assumptions  that  all 
fission  neutrons  are  bom  in  the  fast  group  and  that  there  is  no  up  scatter  from  the  thermal 
group.  The  final  form  of  the  two  energy  group  diffusion  equation  becomes 


V»DiV^j  +  "XiKi  ,  (^1  ^/1  ^1  ^^2  X J2  ^2  ) 


^2  '^a2  ^2  ~  2s12  ^1 


(8) 

(9) 


where 

subcripts  1  and  2  refer  to  groups  1  and  2  respectively 

~  ^rl  ~^S12 

Esi2  =  Cross  section  for  scatter  from  group  1  to  group  2. 
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Two  Dimensional,  One  Energy  Group  Model 

To  develop  this  model,  I  chose  to  use  a  criticality  search  technique  outlined  in  several 
references  (Duderstadt  and  Hamilton,  1976:214-226),  (Clark  and  Hansen,  1964:175-178), 
(Glasstone  and  Sesonske,  1981:208-213),  and  (Ott,  1989:349-356).  This  section  will 
derive  the  methodology  for  the  criticality  search  and  the  numerical  development  to  solve 
the  two  dimensional,  one  energy  group  diffusion  equation. 

As  stated  earlier,  the  energy  level  of  neutrons  within  a  typical  PWR  ranges  from 
about  0.001  eV  to  10  MeV.  Modeling  the  reactor  proves  very  complicated  when 
attempting  to  incorporate  the  entire  energy  range.  Historically,  attempts  to  solve  the 
diffusion  equation  assumed  all  neutrons  were  at  the  same  energy  level.  The  key  to 
solving  the  one  group  model  is  selecting  the  appropriate  macroscopic  cross  section  data. 
The  cross  sections  are  dependent  upon  the  neutron  energy  level.  By  choosing  the 
appropriate  cross  section,  the  one  group  model  can  provide  quantitative  as  well  as 
qualitative  analysis  of  the  reactor  behavior.  The  parameters  chosen  for  this  model  were 
based  on  the  homogenized  data  from  a  typical  reactor  (Duderstadt  and  Hamilton, 
1976:210).  Certainly,  nuclear  reactor  designers  would  not  use  the  one  group  model  for 
design  purposes.  The  value  of  using  a  one  group  model  is  its  ease  of  calculation  and 
proven  qualitative  similarity  to  more  rigorous  models. 

Criticality  Solution  Technique. 

Determining  the  flux  at  criticality  becomes  an  eigenvalue  problem  as 

-V  •  DV(^(r) +2>(r)  =  |v2^<z>(r)  (10) 

rC 
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where  —  is  the  eigenvalue.  For  criticality,  we  seek  k  equal  to  one.  Rewritten  in  matrix 
k 

form,  equation  (10)  yields 

M^  =  -F^  (11) 

k 


where 

F  =  vLj-  (1/cm) 

and  the  operator  M  =  -V  •  DV  +  (1/cm). 

To  solve  this  problem,  we  guess  an  initial  "source"  term  S(r)  and  k  value  where 

Sir)  =  F(/>ir)  =  (r)  and  k  ^  (12) 

and  solve  for  the  flux  using  a  tridiagonal  solver. 

k 

After  solving  for  the  flux  we  must  recalculate  the  source  and  k  values.  The  source  is  easy 
to  recalculate  based  on  known  values. 

(14) 

The  iterative  scheme  is  shown  in  equations  (15)  and  (16).  This  repetitive  process  yields 
the  flux  at  successive  values  until  equation  (17)  is  approximately  true  within  set 
tolerances. 


(15) 

(16) 

(17) 

To  solve  for  the  next  k  iterative  value,  we  recognize 
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(18) 


Mr'-”=^Fr 


^(»1) 


Solving  for  ,  we  then  integrate  the  flux  over  space.  This  is  essentially  averaging  the 
values  to  obtain  a  new  eigenvalue,  where 


An+l) 


The  integration  is  accomplished  numerically  using  the  composite  trapezoid  rule 


J'5(rMr=^  Sia)  +  S{b)  +  2'2sir,) 

where  n  is  the  number  of  mesh  points.  The  iteration  process  continues  until  the 
tolerances  for  k  and  the  source  are  within  a  specified  tolerance. 


An) 


(fj  and/or 


g(n) 


:in) 


{£2 


(19) 


(20) 


(21) 


where 

£1  =0.00001 

£2  =0.015 

The  tolerance  setting  for  is  critical  to  achieving  low  relative  errors  compared  to  the 
mathematical  solution.  Ott  recommends  a  tolerance  of  lE-5  for  most  calculations  (Ott, 
1986;  351).  See  Chapter  III,  Program  Validation  for  details.  As  the  number  of  iterations 
gets  large,  we  expect  the  flux  to  converge  to  the  fundamental  eigenfunction  (Duderstadt 
and  Hamilton,  1976:216-219).  This  will  provide  the  correct  flux  mode  shape  to  enable 
power  and  flux  calculations.  Figure  4  is  a  flowchart  of  the  technique  used  in  the  code  to 
solve  for  the  flux  and  criticality  based  on  the  core  material  composition  and  geometry 
input. 
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Figure  4  Criticality  Search  Technique,  2D  Model 


One  must  provide  an  initial  guess  for  k  and  the  source  flux  in  order  for  the  power 
iteration  process  to  converge  resulting  in  criticality.  Because  it  is  difficult  to  guess  a 
sufficiently  close  guess,  one  must  use  perturbation  theory.  Perturbation  methodology 
assists  the  user  by  adjusting  the  macroscopic  cross  section  until  criticality  is  met. 
Changing  the  macroscopic  cross  section  by  some  small  amount  such  as 

5:,(r)  =  2:,(r)  +  ^,(r),  (22) 

where 

is  the  value  perturbed  by  some  small  positive  or  negative  change  SL^ , 
yields  a  revised  equation  in  matrix  form. 
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M'(^'  =  KF<f 

Jc 


(23) 


The  perturbation  in  the  cross  section  changes  the  diffusion  operator. 

M’=M+SM’  (24) 

where 


^M’=<E„(r). 


We  then  calculate  the  change  in  k  by  applying  the  scalar  product  equation  (23)  with  the 
adjoint  flux  f  of  the  unperturbed  core  obtaining  equation  (25). 

^’) +  {</>*, SM^')  =  j{f,F^’)  (25) 


Using  the  inner  product  of  the  adjoint  operator,  yields 

F*(t>\0^  = 

where,  for  the  one  group  diffusion  model 

^*,F*,  andM*  are  the  adjoint  values 

F=F* 

M=M*. 


(26) 


Substituting  equation  (26)  into  equation  (25)  yields 


^j__r 

^k’  k 


{f,F^')  ■ 


(27) 


However,  this  requires  us  to  know  the  adjoint  and  perturbed  fluxes  that  cannot  be 
calculated  directly.  We  can  rewrite  the  left-hand  side  in  terms  of  reactivity. 


Ap  = 


1_J_ 

k  k 


{f,F<!>’) 


(28) 
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Using  perturbation  theory  we  can  translate  the  unknowns  into  known  values.  A  small 
change  in  absorption  cross  section  is  assumed  to  result  in  a  small  change  in  flux 
(Duderstadt  and  Hamilton,  1976:223).  Expanding  equation  (28)  provides  the  following. 


Ap  = 


{$•  .SLM) 


+  ... 


(29) 


Using  the  self  adjointness  of  the  flux  provides  the  form  required  to  calculate  the  small 
change  in  macroscopic  cross  section, 

U(rfSL^ir)d^r  SZ 

Ap  =  -4 - ; - —  =  — - 

|<?>(r)VS^(r)jV 


(30) 


Vl, 


(31) 


where,  for  criticality, 


k  =1. 


With  a  known  change  in  the  macroscopic  absorption  cross  section,  the  program  user  can 
iterate  the  program  until  criticality  is  achieved  for  the  geometry  and  material  composition 
specified,  if  achievable.  The  change  in  the  macroscopic  absorption  cross  section  can  only 
be  accomplished  physically  by  changing  the  material  composition  in  the  homogeneous 
reactor  core  because 


02) 

1=1 

where 

n  =  the  number  of  materials 

Nf  -  the  number  density  of  the  material  i  (neutrons/cm^ ) 

(7^  =  the  microscopic  absorption  cross  section  of 
material  i  (cm^). 
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Numerical  Development. 

The  numerical  development  of  the  one  group  model  is  based  on  the  central  difference 
approximation  to  the  diffusion  equation  in  right  circular  cylinder  coordinates. 


(I>\r)  =  --(l>\r)+^  (t>{r)  -  ^  v  <zi(r) 

r  D  D 


(33) 


Expanding  </>  in  a  Taylor  series  about  r 


dr 


2  j2. 


+  - 


A"  d 


2  dr^ 


(34) 


dr 


2  dr^ 


(35) 


Now  adding  equations  (34)  and  (35)  yields  the  standard  central  difference  formula  with 
an  order  of  A^ . 


d^(t> 


dr^ 


For  the  standard  differential  the  central  difference  yields 


dr 


_  Ym  Vi-\ 

2A 


The  final  form  of  the  numerical  equation  becomes 


h\ 


1-P(0-  (|>i.l+{2  +  h^q{ri)<p^+  -l  +  p(r)- 


(l>M  =  -h^r(r) 


(36) 


(37) 


(38) 


where 
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h  is  the  distance  between  nodes  (cm) 

p(r)=-  (1/cm) 

i; 


q(r)=^  (1/cm^) 

, ,  1  ^  ,  1  neutrons -cm. 

^ - )• 

kD  cm  cm  -sec 


The  boundary  conditions  are  shown  in  Figure  5  for  a  typical  reactor  core. 


Figure  5  Boundary  Conditions 

Using  a  standard  tridiagonal  solver  rapidly  provides  the  flux  values  for  this  system  of 
equations  along  the  radius  of  the  core. 
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Three  Dimensional,  One  Energy  Group  Model 

I  used  the  same  criticality  search  technique  for  this  model  as  in  the  two  dimensional 
model;  however,  the  numerical  solution  technique  is  quite  different.  Adding  the  third 
dimension  increased  the  boundary  conditions  and  allowed  me  to  develop  a  model  that 
allows  the  user  to  choose  to  calculate  the  power  distribution  and  criticality  in  either  a  half 
or  a  quarter  of  the  reactor  core.  Duderstadt  recommends  using  a  Gauss-Seidel  or  a 
successive  relaxation  method  to  solve  the  numerical  equations  (Duderstadt  and  Hamilton, 
1976:191).  I  chose  to  develop  and  use  a  blocked  tridiagonal  solver  because  it  reduced  the 
computational  time  and  computer  memory  requirements  over  those  recommended. 

Criticality  Solution  Technique 

The  three  dimensional  model  uses  the  same  criticality  iterative  search  technique  to 
solve  for  the  flux  as  the  two  dimensional  model;  however,  the  derivation  of  the  volume 
source  integration  used  in  equation  (19)  to  solve  for  is  more  complicated.  See 
Appendix  A  for  a  complete  derivation.  As  shown  in  Figure  6,  the  three  dimensional 
model  uses  the  blocked  tridiagonal  solver  vice  the  tridiagonal  solver  for  the  two 
dimensional  model. 
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Figure  6  Criticality  Solution  Technique 


Numerical  Development. 

In  three  dimensions,  the  diffusion  equation  in  right  circular  cylinder  coordinates 
becomes 


rdr^  dz\ 


D 


D 


where  the  appropriate  Laplacian  is 

= - r — H - + - . 

r  dr  dr  r^  d&^  dz^ 


(39) 


Owing  to  symmetry 


=  0. 


Using  central  differences  and  collecting  terms  yields 
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(40) 


Ar" 


r2Ar 


+ 


-2^,. ; 


and 


1  1 


Ar^  2rAr 


-2  2  s. 


^Ar'  Az' 


Az" 


D  D 


1  ^ 


Ar^  2rAr  ^ 


<^w,y 


(41) 


where  the  boundary  conditions  are  shown  in  Figure  7.  The  model  provides  solutions  for 
half  of  the  reactor  core  or  quarter  of  the  reactor  core  owing  to  symmetry  of  the 
homogeneous  system. 


core  height 
+Z  P(i,m)=0 


-Z 

Boundary  Conditions  Upper  Quarter  Rx 


core  height 
+Z  P(i,m)=0 


P(i,0)=0 


Boundary  Conditions  Right  Half  Rx 


Figure  7  Boundary  Conditions 


To  explain  the  method  of  converting  the  numerical  equation  into  a  system  of  linear 
equations,  I  will  use  the  sample  mesh  system  shown  in  Figure  8. 
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Figure  8  Example  4x4-Three  Dimensional  Mesh  System 

To  convert  these  equations  into  a  solvable  linear  algebra  system,  it  is  necessary  to 
convert  the  (i,j)  indices  into  a  single  value  using 

l  =  i  +  im-l-j)in-V)  (42) 

where  m  is  the  number  of  nodes  along  the  z-axis  and  n  is  the  number  of  nodes  along  the 
radius  (Barden  and  Faires,  1997:676).  This  conversion  results  in  re-numbering  the 
interior  mesh  points  as  shown  in  Figure  9. 
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Figure  9  Relabeled  Interior  Mesh  Points 

Referring  to  equation  (41),  this  system  has  five  terms  instead  of  the  three  terms  used 
in  the  two-dimension  model.  This  corresponds  to  the  sample  9x9  matrix  shown  in 
equation  (43).  The  matrix  is  a  tridiagonal  system  with  a  sub  and  super  diagonal.  The  sub 
and  super  diagonal,  in  this  case,  contain  constant  and  equal  values  in  each  component. 
The  upper  and  lower  tridiagonal  diagonals  are  variables  that  depend  on  the  position  along 
the  radius  as  shown  in  equation  (41). 
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(43) 


This  same  pattern  applies  to  any  size  matrix  depending  upon  the  number  of  interior 
mesh  points  selected.  The  boundary  conditions  are  incorporated  into  positions  Pi,  P4,  and 
P7  in  the  matrix. 

Matrix  Solution  Methods. 

There  are  several  methods  available  to  solve  these  types  of  matrix  problems.  One  is 
the  Jacobi  method.  This  method  converges  too  slowly  for  practical  use  in  large  matrices 
because  of  the  number  of  required  operations.  Another  method,  Gauss-Seidel,  is  often 
used  to  solve  small  to  medium  sized  matrices.  Gauss-Seidel  also  converges  slowly  and, 
like  the  Jacobi  method,  requires  storing  every  point  in  the  matrix  in  computer  memory. 
As  a  result,  it  is  slow  and  computationally  inefficient.  Successive  over-relaxation  (SOR) 
is  an  improved  version  of  the  Gauss-Seidel  method.  It  makes  an  over  correction  by 
anticipating  future  corrections.  To  reduce  the  error  by  a  factor  of  10'^,  the  SOR  method 
requires  on  the  order  of  J  iterations  compared  to  J^  for  the  other  methods  (Press, 
1996:858)  where 


and 


r^-pj 

r  =  the  rth  stage  of  the  iteration  process 
J  =  the  number  of  iterations. 


(44) 
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Each  of  these  methods  requires  filling  and  using  the  entire  matrix  to  obtain  a  solution. 
This  is  computationally  inefficient  when  dealing  with  relatively  large  matrices.  Using  a 
reactor  size  of  120  cm  radius  and  360  cm  height  and  0.5  cm  node  spacing  requires  a 
171,841  X  171,841  matrix.  This  is  a  fairly  large  matrix  and  using  any  of  the  above 
techniques  would  increase  the  computation  time  and  memory  requirements  of  the 
computer. 

To  solve  this  system,  I  chose  a  blocked  tridiagonal  solver  technique.  This  method 
requires  only  storing  the  values  in  the  diagonals  of  the  tridiagonal  matrix  of  size  (m-1  x 
n-1)  as  compared  to  (m-1  x  n-1)^  for  the  Gauss-Seidel  method.  Using  the  previous  9x9 
example,  the  system  of  equations  in  (43)  becomes 


where  D  is  the  sub  and  super  diagonal  with  constant  coefficients,  Ai,2,3  are  tridiagonal 
matrices,  Xi,2,3  are  the  unknown  flux  values  at  the  interior  mesh  points,  and  61,2,3  are  the 
solutions  at  each  interior  mesh  point. 

This  can  be  further  broken  down  into  three  equations  that  can  each  be  solved  using  a 
standard  tridiagonal  solver. 

DX,  -I- \X^  +  DX3  =  (46) 

DX,  +  A,X,=B, 
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Inverting  the  D  diagonal  (which  is  constant)  and  rearranging  the  equations  allows 
them  to  become 

D-^Aj  Xf  ’  =  D-^Bj  -  X”!j  -  X (47) 

where)  is  the  jth  row  along  the  z-axis.  This  allows  one  to  solve  the  system  only  using  the 
main  tridiagonal  components.  To  solve  the  system  of  equations  in  (47),  set  the  initial 
values  of  Xj  to  zero  and  solve  each  equation  in  a  tridiagonal  solver  making  use  of  the 

previous  X^.  value.  This  iterative  approach  is  similar  to  the  Gauss-Seidel  approach 
without  the  excess  storage  or  computations.  Once  the  values  of  Xj  are  within  the  set 
tolerance  of  the  previous  X^.j ,  then  the  process  has  converged  to  the  solution.  This 

tolerance  level  is  critical  to  achieving  low  relative  errors  between  the  old  and  the  new 
flux  values.  For  example,  setting  the  tolerance  equal  to  0.001  provides  a  maximum 
relative  error  of  18  %  for  the  axial  power  using  a  mesh  spacing  of  one  centimeter. 
Changing  the  tolerance  to  lE-6,  reduces  the  maximum  relative  error  to  8  %.  Reducing 
the  tolerance  does  however  increase  the  computational  time  dramatically.  See  Chapter 
III,  Program  Validation  for  further  details. 

Three  Dimensional,  Two  Energy  Group  Model 

Criticality  Solution  Technique 

The  solution  technique  is  the  same  as  the  one  group  homogeneous  method  except  one 
must  guess  an  initial  value  for  <p^  and  ^2  to  solve  for  in  equation  (8).  is  then  used  to 
solve  for  (p^  in  equation  (9).  The  code  iterates  as  in  the  one  group  model  solving  for 
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and  per  equations  (16)  and  (19)  and  then  checks  for  convergence  per  equation 
(21).  This  is  shown  in  Figure  10. 


Figure  10  Criticality  Solution  Technique  for  Two  Energy  Groups 


There  is  a  difference  in  calculating  the  perturbation  of  the  cross  sections.  The  method 
is  not  as  simple  as  in  the  one  energy  case  because  the  multi-group  criticality  problem  is 
not  self-adjoint.  I  did  not  include  the  perturbation  because  of  the  complexity  of  having  to 
change  multiple  parameters  in  both  energy  groups.  The  method  is  outlined  in  Duderstadt 
and  Hamilton. 

Numerical  Development. 

The  development  of  the  numerical  solution  technique  is  very  similar  to  the  three 
dimensional  one  group  technique;  however,  there  are  now  two  coupled  equations  to 
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solve.  The  first  equation,  derived  from  applying  central  differences  to  equation  (8),  is 
shown  in  equation  (48). 


1 
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where 

■S'  =  (vj  X  fx  +  ^2  X  f2  ^2  ) 

Z.=l- 

The  program  solves  this  equation  for  the  fast  flux,  ,  and  then  solves  the  second  coupled 
equation  for  the  thermal  flux,  ^2  using  the  blocked  tridiagonal  solver.  The  second 
coupled  equation,  derived  by  applying  central  differences  to  equation  (9),  is 
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Ifl.  Program  Development  and  Validation 
Program  Development 

In  general,  the  program  provides  much  flexibility  for  the  user.  The  user  can  choose 
between  a  two  dimensional,  one  energy  group  model,  a  three  dimensional  one  energy 
group  model,  or  a  three  dimensional  two  energy  group  model.  The  program  provides  a 
drop  down  Windows  menu  to  allow  the  user  to  select  various  command  options  such  as 
models  to  run  or  help  files  to  view.  Within  each  model,  the  user  can  point  the  mouse 
over  input  boxes  and  get  definitions  or  additional  explanations.  This  serves  to  assist  the 
user  in  understanding  either  the  physics  involved  or  how  to  proceed. 

The  program  requires  initial  data  input  that  represents  a  homogenized  reactor  core  in 
cylindrical  geometry.  It  requires  the  user  to  input  the  reactor  core  dimensions  and 
material  composition/cross  sections.  If  the  user  is  not  sure  of  the  material  composition, 
the  program  provides  recommended  input  values.  The  values  for  the  macroscopic  cross 
sections  (absorption  and  fission),  the  diffusion  coefficient,  the  number  of  interior  mesh 
points,  and  the  initial  guess  values  for  Ineffective  and  flux  are  required  for  input.  The 
program  will  only  allow  for  equal  mesh  spacing  in  both  the  axial  and  radial  directions 
and  requires  the  core  height  and  radius  be  multiples  of  mesh  spacing.  The  boundary 
conditions  are  established  within  the  program.  For  all  exterior  points  the  flux  is  assumed 
to  be  zero.  At  the  center  of  the  reactor  core,  the  flux  is  set  equal  to  the  flux  value  at  the 
first  interior  mesh  point  away  from  the  centerline.  This  makes  use  of  the  symmetry  of  the 
homogeneous  core  and  meets  the  requirement  that  the  first  derivative  of  the  flux  equals 
zero  at  the  center. 
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Given  the  initial  input,  the  program  will  calculate  the  initial  source  term 
value  (vS/  for  the  right  hand  side  of  the  diffusion  equation  at  each  interior  mesh  point. 

If  the  user  requested  to  solve  the  two  dimensional  simulation,  the  program  calls  a  linear 
finite  difference  solver  that  uses  the  Grout  factorization  method  for  tridiagonal  linear 
systems.  The  solver  returns  a  vector  of  flux  values  along  the  radial  interior  mesh  points. 
Next,  the  program  begins  to  iterate  to  converge  to  the  solution.  It  evaluates  the  tolerance 
between  the  old  and  the  new  flux  values  and  the  old  and  new  keffective  values.  To  calculate 
the  new  keffective>  the  program  integrates  the  old  and  new  source  terms  over  space  and 
essentially  averages  them  as  in  equation  (19).  With  the  updated  source  and  keffective,  the 
program  iterates  again.  This  process  continues  until  convergence  is  met.  Next,  the 
program  checks  if  the  reactor  is  critical.  If  keffective  is  greater  than  or  equal  to  one,  the 
system  is  critical.  If  the  reactor  is  not  critical,  the  program  uses  perturbation  theory  to 
provide  a  revised  macroscopic  cross  section  that  should  assist  the  user  in  achieving  a 
critical  reactor.  In  either  case,  the  data  are  automatically  loaded  onto  an  Excel 
spreadsheet  and  plotted  in  Excel.  The  Excel  chart  is  automatically  updated  onto  the 
Visual  BASIC  form  and  saved  to  a  location  provided  by  the  user.  The  user  of  the 
program  does  not  see  Excel  running  in  the  background. 

For  the  three  dimensional  one  energy  group  problem,  the  user  has  the  option  of 
selecting  to  calculate  the  power  profile  for  either  a  half  or  a  quarter  of  the  reactor  core. 
Because  of  the  symmetry  of  the  core  geometry,  the  solution  to  the  half  of  the  reactor  core 
is  a  mirror  image  of  the  quarter  core  solution.  The  iterative  processes  are  very  similar  to 
the  two  dimensional  problem;  however,  the  solver  is  different.  Because  the  numerical 
analysis  problem  generates  a  blocked  tridiagonal  system,  the  solution  technique  changes. 
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The  program  builds  blocked  tridiagonal  vectors  that  create  a  (m-1)  system  of  equations. 

It  fills  each  diagonal  component  with  its  corresponding  coefficient  and  then  iterates 
through  the  Grout  factorization  method  for  each  tridiagonal  linear  system  until  the  desired 
convergence  is  met  for  the  iterative  solution.  Next,  the  program  updates  keffective  and  tests 
for  convergence  of  the  old  and  new  source  terms  similar  to  the  two  dimensional  problem. 
If  the  reactor  is  not  critical,  the  program  will  again  recommend  an  adjusted  value  for  the 
macroscopic  cross  section  and  plot  both  the  radial  and  axial  power  values  in  Excel. 

The  three  dimensional,  two  energy  group  model  uses  the  same  process  as  the  three 
dimensional,  one  energy  group  model  with  a  few  modifications.  In  addition  to  the 
previous  input  requirements,  the  user  must  provide  the  macroscopic  scatter  and  removal 
cross  sections  for  the  fast  and  thermal  energy  ranges  as  well  as  the  fast  and  thermal  flux 
initial  guesses.  Using  this  data,  the  program  solves  for  the  fast  flux  using  the  blocked 
tridiagonal  solver  and  then  it  uses  that  value  to  solve  for  the  thermal  flux,  again  using  the 
blocked  tridiagonal  solver.  The  program  iterates  as  before  until  convergence  is  met  for 
the  total  source  term  values  and  k. 

The  user  can  double  click  each  Excel  chart  on  the  form  and  open  Excel  to  access  the 
chart  or  the  data.  A  complete  program  along  with  typical  output  data  is  included  in 
Appendices  C,  D,  E,  and  F. 

Operating  the  Code 

Each  window  in  the  program  is  a  form  that  allows  the  program  user  several  options. 
There  is  typically  a  drop  down  menu  window  structured  as  shown  in  Figure  11. 
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Figure  11  Typical  Drop  Down  File  Menu 


The  initial  welcome  form  is  shown  in  Figure  12  below. 
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Figure  12  Welcome  Screen 


The  PrintForm  menu  function  prints  the  current  form  to  the  default  printer.  The  Exit 
menu  function  exits  the  program.  Clicking  the  View  menu  function  provides  two 
options.  The  program  user  can  either  view  a  schematic  of  the  reactor  or  the  boundary 
conditions  as  shown  in  Figure  13  and  Figure  14  respectively.  From  either  of  them  the 
user  can  return  to  the  previous  screen  by  clicking  the  return  menu.  To  choose  one  of  the 
three  models  to  run,  click  the  Select  Model  menu  and  choose  the  desired  model.  From 
the  chosen  model,  the  user  can  return  to  the  starting  window  or  select  from  any  of  the 
options  shown.  Clicking  the  run  menu  function  will  execute  the  selected  module.  The 
help  menu  function  will  bring  up  a  window  with  a  help  object  written  in  Word.  To  view 
the  help  information,  double  click  the  object  window  and  scroll  through  the  file.  To  close 
the  file  and  return  to  the  previous  screen,  click  the  “x”  on  the  windows  screen. 
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Boundary  Conditions 


Figure  14  Boundary  Condition  Form 
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The  program  has  several  built  in  error  commands  that  prevent  the  program  from 
crashing  should  the  user  fail  to  input  required  data  or  input  incorrect  data.  For  example, 
if  the  user  fails  to  input  ,  the  program  will  anticipate  and  prevent  the  division  by 

zero  in  the  code.  The  code  will  display  a  message  explaining  the  error  and  insert 
^effecihe  =  b.9  .  Whilc  thc  codc  is  certainly  not  yet  totally  failsafe,  it  provides  error 

corrections  to  many  anticipated  runtime  errors. 

To  run  the  two  dimensional  one  energy  group  model,  the  code  requires  input  for  2^, 
v2f,  diffusion  coefficient,  radius,  and  All  appropriate  data  must  include  units 
of  centimeters.  Each  of  the  three  reactor  model  forms  provides  tips  for  the  user  when  the 
cursor  is  placed  over  some  of  the  input  description  boxes.  The  user  must  also  input  the 
mesh  spacing  between  mesh  points.  For  typical  reactors  a  mesh  spacing  of  0.5  cm 
provides  a  maximum  relative  error  of  less  than  two  percent  and  runs  within  a  few 
seconds.  Finally,  the  user  must  input  a  file  name  and  location  to  save  the  Excel 
workbook  with  an  extension  of  filename.xls. 

Figure  15  shows  an  example  of  the  two  dimensional  model  form.  The  form  lets  the 
user  know  if  the  reactor  is  critical  and  if  not  provides  a  recommended  2^  that  will  drive 
the  reactor  to  criticality. 
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Figure  15  2D,  One  Energy  Group  Form 


It  also  provides  a  power  distribution  plot.  To  access  the  plot  and  data,  double  click  the 
chart  and  Excel  will  open  the  workbook  containing  the  data  and  chart.  Closing  Excel 
returns  the  user  to  the  form;  however,  the  chart  will  not  be  resized  to  fit  the  screen.  The 
chart  will  automatically  resize  upon  running  another  problem. 

To  run  the  three  dimensional  one  energy  group  model,  input  the  same  information 
required  for  the  two  dimensional  model  along  with  the  reactor  height.  The  program 
requires  the  user  to  select  whether  to  calculate  the  power  distribution  for  either  a  half  or  a 
quarter  of  the  core.  The  mesh  spacing  will  be  equally  applied  to  the  axial  and  radial 
directions.  Additionally,  there  is  an  option  to  choose  between  reduced  relative  error; 
slower  run  times  and  average  relative  error;  faster  run  times.  These  correspond  to  the 
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results  provided  in  Table  5,  Table  6,  and  Table  10.  Figure  16  shows  the  results  of 

running  a  sample  problem  using  half  of  the  core. 

1^— .aiiipiai 


Figure  16  3D,  One  Energy  Group  Form 


The  reactor  core  height  input  corresponds  to  the  height  of  either  the  half  or  quarter 
reactor  core  selection.  Again,  the  program  will  provide  the  criticality  information. 

The  three  dimensional  two  energy  group  model  requires  more  information  than  the 
previous  models.  The  input  data  for  both  the  thermal  and  fast  energy  groups  must  also 
include  ^scatter.  •  Table  2  provides  sample  input  data  for  two  energy  groups. 

Additionally,  the  user  must  select  the  radial  and  axial  positions  to  plot  the  power 
distribution.  Figure  17  displays  the  input  box  to  select  the  radial  position  upon  which  to 
plot  the  axial  power  distribution.  The  maximum  axial  power  will  occur  at  the  zero 
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position  for  either  the  half  or  quarter  core  selection.  For  the  radial  power  distribution,  the 
maximum  will  occur  at  the  zero  axial  position  for  the  quarter  core  or  center  axial  position 
for  the  half  core.  The  plots  are  independent  of  one  another  allowing  the  user  to  select  the 
power  distribution  along  any  section  of  the  core. 

Like  the  three  dimensional  one  energy  group  model,  the  user  must  choose  between 
reduced  relative  error;  slower  run  times  and  average  relative  error;  faster  run  times. 


Figure  17  3D,  Two  Energy  Groups  Axial  Choice 


Figure  18  is  an  example  of  the  final  solution.  The  plots  include  the  thermal,  fast,  and 
total  power  distributions.  As  in  the  other  models,  the  program  user  can  access  the 
Excel  workbook  directly  by  double  clicking  either  chart.  The  workbook  will  contain 
the  power  distribution  data  for  the  radial  and  axial  plots  on  separate  worksheets. 
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Figure  18  3D,  Two  Energy  Group  Form 
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Program  Validation 

To  validate  the  code,  I  tested  each  of  the  three  models  independently.  I  compared  the 
code’s  normalized  power  distribution  to  the  normalized  analytical  solution  to  determine 
the  point-by-point  relative  error.  For  the  three  dimensional  models,  I  tested  and  validated 
the  blocked  tridiagonal  solver  using  various  matrices  solved  with  Mathamatica  before 
testing  the  entire  code. 


Two  Dimensional,  One  Energy  Group  Model. 

I  tested  the  model  using  data  for  a  typical  homogeneous  reactor  as  shown  in  Table  1. 


Table  1  One-Speed  Reactor  Input  Data 


Reactor  Data 

0.1532  1/cm 

0.1570  1/cm 

Diffusion  Coefficient 

9.21  cm 

0.0362  1/cm 

I  compared  the  model  to  the  flux  distribution  profile  for  an  infinite  right  circular  cylinder. 


(p  —  Jq 


(50) 


where 


Vq  -  2.405  =  smallest  zero  of  Jq 
R  is  the  radius  of  the  cylinder 
r  is  the  position  along  the  radius. 
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Figure  19  shows  the  normalized  flux  profile  for  the  reactor  with  a  radius  of  120  cm. 
Figure  20  shows  the  normalized  flux  profile  for  the  data  from  the  one-group  model. 
Figure  21  is  a  combination  of  both.  They  overlap  each  other  indicating  that  the  one- 
speed  model  is  producing  the  correct  fundamental  flux  mode  shape.  Using  a  radius  of 
120  cm  and  a  mesh  spacing  of  0.5  cm,  the  maximum  relative  point-by-point  error  was 
0.02. 


Figure  19  Flux  profile  of  infinite  right  circular  cylinder  (Bessel  J  function) 
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Normalized  Flux  I  Normalized  Flux 


Figure  20  Normalized  flux  from  one-speed  computer  model 


Two  Dimensional,  One  Energy  Group 
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Figure  21  Flux  proflie  of  Bessel  J  function  and  one-speed  computer  model  together 
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To  test  the  two  models,  I  verified  the  blocked  tri diagonal  solver  and  then  each  code 


separately.  I  initially  tested  the  blocked  tridiagonal  solver  using  Mathematica  with  a 
diagonally  dominant  system  of  equations  as  shown  in  the  augmented  matrix  (51). 

■-5.16  1.25  0  1  0  0  0  0  0  -1.894‘ 

0.5  -5.16  1.167  0  1  0  0  0  0  -1.894 

0  0.75  -5.16  0  0  1  0  0  0  -1.894 

1  0  0  -5.16  1.25  0  1  0  0  -1.894 

0  1  0  0.5  -5.16  1.167  0  1  0  -1.894  (51) 

0  0  1  0  0.75  -5.16  0  0  1  -1.894 

0  0  0  1  0  0  -5.16  1.25  0  -1.894 

0  0  0  0  1  0  0.5  -5.16  1.167  -1.894 

0  0  0  0  0  1  0  0.75  -5.16  -1.894 

The  approximate  solution  is 

100000000  0.67868f 
010000000  0.678681 
00  1  0  0  0  0  0  0  0.678681 
0  0  0  1  0  0  0  0  0  0.678681 

000010000  0.678681  .  (52) 

000001000  0.678681 
000000100  0.678681 
0  0  0  0  0  0  0  1  0  0.678681 

0  0  0  0  0  0  0  0  1  0.678681 

In  this  case,  the  system  resulted  in  a  nine  by  nine  matrix  where  m=4  and  n=4.  The 
pattern  and  coefficients  are  similar  to  the  actual  pattern  and  coefficients  generated  by 
numerically  solving  the  diffusion  equation  using  the  one  group  data.  The  solution 
generated  by  Mathematica  matched  the  solution  given  by  the  blocked  tridiagonal  solver. 
The  blocked  tiidiagonal  solver  typically  converges  to  the  approximate  solution  within  10- 
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15  iterations  with  a  tolerance  of  lE-6,  using  mesh  spaces  greater  than  10  centimeters. 
Their  agreement  indicates  that  the  solver  is  converging  to  the  correct  solution. 

Varying  the  values,  constant  and  non-constant  coefficients,  of  the  lower  and  upper 
diagonals  while  maintaining  diagonal  dominance  resulted  in  correct  solutions  as  well. 
This  was  one  of  the  several  comparisons  I  made,  each  of  which  agreed  with  the  blocked 
tridiagonal  solver’s  solution. 

Next,  I  verified  that  the  codes  produced  the  approximate  correct  power  distribution 
within  the  reactor  core.  The  compiled  program  was  tested  using  an  AMD-K3  450 
megahertz  personal  computer  with  64  megabytes  of  random  access  memory.  I  compared 
the  relative  maximum  error  of  the  normalized  numerical  solutions  at  each  mesh  point  to 
the  normalized  mathematical  solution  for  a  finite  cylinder. 


Hr,z)  =  Jc 


2.405r  (  nz 

-  Cos  — 

i?  \  n 


where 

H  =  the  height  of  the  cylinder 
z  =  the  position  along  the  z  axis  (height) 

For  the  three  dimensional  one  group  model,  I  used  the  same  input  data  used  in  the 
two  dimensional  one  group  model.  I  used  the  data  in  Table  2  as  input  for  the  three 
dimensional  two  energy  group  model  (Duderstadt  and  Hamilton,  1976:312). 


45 


Table  2  Two  Group  Diffusion  Theory  Constants  (Table  7-2  Duderstadt) 


Group 

Constant 

Group  1  Fast 

Group  2 
Thermal 

Ea  (1/cm) 

0.01207 

0.1210 

vEf  (1/cm) 

0.008476 

0.18514 

Ef  (1/cm) 

0.003320 

0.07537 

Diffusion 

Coefficeint 

(cm) 

1.2627 

0.3543 

SigmaRemoval 

(1/cm) 

0.02619 

0.1210 

SigmaScatter 

(1/cm) 

0.01412 

0 

I  originally  limited  the  blocked  tridiagonal  solver  to  20  iterations  because  of  the 
typical  convergence  within  10-15  iterations.  However,  this  produced  erroneous  results 
as  the  mesh  spacing  was  reduced  to  less  than  six  centimeters.  The  blocked  tridiagonal 
solver  failed  to  completely  converge  after  20  iterations  causing  the  large  relative  errors 
shown  in  Table  3  and  Table  4.  Notice  the  relative  error  converges  and  then  begins  to 
diverge  below  mesh  spacings  of  around  4-6  centimeters.  Table  3  and  Table  4  show  the 
relative  errors  and  approximate  code  running  times  for  several  test  runs  for  half  and 
quarter  of  the  core  test  runs. 
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Table  3  Relative  Errors  for  Quarter  Core  3D  Models,  Blocked  Tridiagonal 
Tolerance  =lE-8  with  Maximum  of  20  Iterations 


One  Energy  Group 


Two  Energy  Groups 


Mesh 

Spacing 

Centimeters 


Radial 

Maximum 

Relative 


Axial 

Maximum 

Relative 


Running 

Time 

(MimSec) 


Radial 

Maximum 

Relative 


Axial 

Maximum 

Relative 


Running 

Time 

(MimSec) 


15_ 

12 

10 

6_ 

_5_ 

3 


Error 

0.1 

0.07 

0.04 

0.03 

0.02 

0.01 

0.008 

0.006 

0.005 


Error 

0.1 

0.06 

0.05 

0.04 

0.03 

0.01 

0.03 

0.2 

0.7 


00:08 

00:13 

00:21 

00:36 

00:52 

02:20 

03:20 

05:27 

Untimed 


Error 

0.1 

0.07 

0.04 

0.03 

0.02 

0.01 

0.008 

0.006 

0.005 


Error 

0.1 

0.06 

0.05 

0.04 

0.03 

0.02 

0.005 

0.1 

0.5 


00:16 

00:21 

00:32 

00:51 

01:21 

03:52 

05:52 

09:54 

Untimed 


Table  4  Relative  Errors  for  Half  Core  3D  Models,  Blocked  Tridiagonal  Tolerance 

=lE-8  with  Maximum  of  20  Iterations 


One  Energy  Group 

Two  Energy  Groups  (Total  Power) 

Mesh 

Spacing 

Centimeters 

Radial 

Maximum 

Relative 

Error 

Axial 

Maximum 

Relative 

Error 

Running 

Time 

(Min:Sec) 

Radial 

Maximum 

Relative 

Error 

Axial 

Maximum 

Relative 

Error 

Running 

Time 

(Min:Sec) 

30 

0.1 

0.1 

0.0001 

20 

0.07 

00:21 

0.07 

0.0001 

15 

0.04 

0.00009 

00:39 

0.0001 

00:59 

12 

0.00009 

01:12 

0.0001 

01:49 

10 

0.02 

■IKIIIIIIIVM 

01:42 

■ilrJElM 

6 

0.01 

■SSEH 

04:41 

5 

0.008 

06:41 

0.008 

0.01 

4 

0.006 

0.6 

3 

0.005 

Untimed 

0.005 

0.4 

I  initially  thought  the  divergence  was  due  to  instability  of  the  finite  central  difference 
method;  however.  Figure  22  shows  that  the  relative  error  was  not  symmetrical. 
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3D,  One  Energy  Group,  Half  Reactor 


Figure  22  Intitial  Test,  Mesh  Spacing  =  3  cm 

This  suggested  that  the  error  was  not  due  to  instability.  Increasing  the  maximum 
number  of  iterations  in  the  blocked  tridiagonal  solver  from  20  to  1000  allowed  for 
complete  convergence  with  smaller  mesh  sizes  and  corrected  the  error. 

The  tolerances  set  for  the  convergence  of  the  multiplication  factor  k  and  the  blocked 
tridiagonal  solver  are  critical  to  achieving  useful  solutions.  For  most  cases,  a  k  tolerance 
of  lE-5  provides  acceptable  results  (Ott,  1989:351).  Table  5  and  Table  6  provide  a 
summary  of  the  results  using  a  tolerance  of  0.001  for  the  blocked  tridiagonal  solver  and  a 
k  tolerance  of  lE-5.  This  tolerance  setting  provided  maximum  relative  errors  of  less  than 
six  percent  for  a  mesh  spacing  of  two  centimeters  when  analyzing  a  quarter  of  the  reactor. 
Analyzing  half  the  reactor  core  increased  the  maximum  relative  error  and  run  times  as 
expected.  I  chose  this  setting  because  it  yielded  reasonable  results  with  fast  run  times. 


Table  5  Relative  Error  for  Quarter  Core  3D,  Blocked  Tridagonal  Tolerance  =  0.001 


One  Energy  Group 


Two  Energy  Groups 


0.04 


0.03 


0.02 


0.01 


9_ 

1 


0.006 


0.006 


0.01 


Axial 

Maximum 

Relative 

Error 

0.2 


Running 

Time 

(MimSec) 

00:06 


Radial 


Axial 


Maximum  Maximum 


Relative 

Error 


0.1 


6 

00:49 

4 

01:32 

bseeh 

5 

1 

No  test 


Relative 

Error 


0.2 


2 _ 

.1 


0.1 

0.1 


0.1 


.09 

.08 


.04 

.02 


No  test 


Running 
Time 
(Min:  Sec) 


00:09 


00:09 

00:10 


00:11 

00:13 


00:33 


P 


03:47 

14:13 


No  test 


Table  6  Relative  Errors  for  Half  Core  3D,  Blocked  Tridiagonal  Tolerance  =  0.001 


What  impact  does  reducing  the  tolerance  of  the  blocked  tridiagonal  solver  have  on  the 
maximum  relative  error?  Table  7  and  Table  8  show  the  results  of  changing  the  tolerance 
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to  lE-4  and  lE-6  respectively  for  half  of  the  core  using  two  energy  groups.  There  was 
not  a  significant  reduction  in  the  maximum  relative  error  by  changing  the  tolerance  to  lE- 
4.  Reducing  the  tolerance  to  lE-6  did  not  significantly  reduce  the  error  in  the  radial 
direction;  however,  the  error  was  reduced  by  over  ten  percent  in  the  axial  direction.  The 
trade  off  is  doubling  the  run  time  from  approximately  29  minutes  as  shown  in  Table  6  to 
60  minutes  as  shown  in  Table  8. 


Table  7  Half  Core  3D,  Blocked  tridiagonal  Tolerance  =  lE-4 


Two  Energy  Groups 

Mesh  Spacing 
Centimeters 

Radial  Maximum 
Relative  Error 

Axial  Maximum 
Relative  Error 

Running 

Time 
(Min:  Sec) 

2 

0.006 

0.2 

33:42 

Table  8  Half  Core  3D,  Blocked  Tridiagonal  Tolerance  =  lE-6 


Two  Energy  Groups 

Mesh  Spacing 
Centimeters 

Radial  Maximum 
Relative  Error 

Axial  Maximum 
Relative  Error 

Running 

Time 
(Min:  Sec) 

2 

0.006 

0.08 

60:05 

Reducing  the  tolerance  to  lE-8  provided  even  better  results  as  shown  in  Table  9. 
This  tolerance  reduced  the  maximum  axial  error  as  shown  in  Table  6  by  approximately 
19  percent  for  the  mesh  spacing  of  one  centimeter. 


Table  9  Half  Core  3D,  Blocked  Tridiagonal  Tolerance  =  lE-8 


Two  Energy  Groups 

Mesh  Spacing 
Centimeters 

Radial  Maximum 
Relative  Error 

Axial  Maximum 
Relative  Error 

Running 

Time 
(Min:  Sec) 

1 

0.01 

0.1 

Several  Hours 
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Figure  23  is  the  normalized  power  distribution  plot  compared  to  the  normalized 
mathematical  solution  for  the  reduced  blocked  tridiagonal  tolerance  of  lE-8.  Although 
the  maximum  error  is  about  10  percent,  the  error  is  symmetric  about  the  center  of  the 
reactor  core. 

3D,  2  Energy  Groups,  Half  Rx 


■—  Numerical  Solution 

- Mathematical 

Solution 


0  40  80  120  160  200  240  280  320  360 
Height  (cm) 

Figure  23  3D,  Total  2  Energy  Groups,  Blocked  Tridiagonal  Tolerance  =  lE-8  with 

Mesh  Spacing  =  1  cm 

To  reduce  the  maximum  relative  error  even  more,  I  set  the  tolerances  for  convergence 
of  the  multiplication  factor  k  and  the  blocked  tridiagonal  solver  to  lE-7  and  lE-8 
respectively.  Table  10  is  a  summary  of  the  maximum  relative  errors  for  test  runs 
calculating  the  power  distribution  for  a  half  reactor  core,  using  two  energy  groups.  The 
data  indicates  that  the  radial  maximum  relative  error  continues  to  reduce  while  the  axial 
error  remains  constant  at  about  0.001.  As  before,  the  running  times  increase  significantly 
as  the  mesh  spacing  reduces  to  one  centimeter. 
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Table  10  Half  Core  3D,  k  Tolerance  =lE-7,  Blocked  Tridiagonal  Tolerance  =  lE-8 


Two  Energy  Groups 

Mesh  Spacing 
Centimeters 

Radial  Maximum 
Relative  Error 

Axial  Maximum 
Relative  Error 

Running 

Time 
(Min:  Sec) 

30 

0.1 

0.001 

00:10 

20 

0.07 

0.001 

00:15 

15 

0.001 

00:25 

12 

0.03 

0.001 

00:40 

10 

0.02 

0.001 

6 

0.009 

0.001 

5 

0.007 

0.001 

07:55 

4 

0.001 

16:59 

3 

0.001 

Appendix  G  contains  plots  comparing  the  relative  errors  as  indicated  in  Table  5  and 
Table  6  using  a  blocked  tridiagonal  tolerance  of  0.001  and  a  k  tolerance  of  lE-5. 

The  final  code  provides  the  user  with  the  options  as  shown  in  Table  11  below.  This 
provides  the  user  with  the  flexibility  to  choose  between  a  level  of  maximum  relative 
errors  and  run  times  corresponding  to  Table  5,  Table  6,  and  Table  10. 


Table  11  Options  for  Three  Dimensional  Models 


Tolerance  for  k 

Tolerance  for  blocked 
Tridiagonal  solver 

Reduced  relative  error, 
Increased  run  time 

lE-7 

lE-8 

Average  relative  error. 
Faster  run  time 

lE-5 

0.001 
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IV.  Conclusions  and  Recommendations 


Conclusions 

The  model  provides  the  power  distribution  of  a  homogeneous  unreflected  reactor  core 
in  two  or  three  dimensions  using  either  one  or  two  energy  groups  for  a  steady  state 
reactor.  Teachers  can  use  the  program  to  augment  fundamental  nuclear  reactor  courses 
by  providing  students  with  an  additional  resource  to  enhance  learning.  The  program 
allows  the  user  to  modify  the  reactor  dimensions  and/or  core  composition  and  see  the 
impacts  on  the  power  distribution  and  criticality  within  the  reactor  core. 

Reducing  the  allowed  acceptable  tolerance  for  convergence  in  the  blocked  tridiagonal 
solver  and  the  multiplication  factor  will  reduce  the  normalized  maximum  relative  error 
for  the  three  dimensional  models.  The  major  trade  off  is  increasing  the  computational 
time.  Using  lE-5  and  0.01  as  the  convergence  tolerances  for  the  multiplication  factor  and 
blocked  tridiagonal  solver  respectively,  the  model  yields  a  power  distribution  with  a 
maximum  relative  error  of  about  four  percent  for  a  mesh  spacing  of  three  centimeters, 
using  two  energy  groups  for  a  quarter  core  calculation.  Using  lE-7  and  lE-8 
respectively,  the  model  yielded  a  maximum  relative  error  of  about  one  half  of  one  percent 
for  the  half  core  calculations.  To  achieve  such  a  low  relative  error,  the  running  times 
increased  from  about  four  minutes  to  48  minutes.  Because  this  is  a  homogeneous  system, 
one  should  take  advantage  of  symmetry  and  calculate  the  power  distribution  in  a  quarter 
of  the  core. 

The  Visual  BASIC  5.0  program  is  completely  exportable  to  most  Windows  based 
personal  computers.  It  automatically  plots  the  power  distribution,  based  on  the  core 
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dimensions  and  composition  input,  using  Excel  and  links  the  chart  to  the  Visual  BASIC 
5.0  form.  Additionally,  it  calculates  the  multiplication  factor  to  determine  criticality. 

Recommendations 

The  code  is  flexible  enough  to  allow  for  future,  user-friendly  improvements  while  the 
finite  central  difference  method  and  criticality  search  technique  are  the  foundation  for 
more  complex  reactor  codes.  The  model  can  be  the  basis  for  adding  a  heterogeneous 
core  and  other  modules  including  thermal-hydraulics,  control  adjustments,  and  depletion. 

Although,  the  code  has  some  built  in  error  checks,  it  is  not  totally  “crash  proof’. 
Several  additional  error  checks  should  be  added  as  the  program  is  used  and  tested  by 
teachers  and  students  alike.  Additionally,  an  improved  help  file  and  automated  read  input 
statement  should  be  added. 

One  approach  to  developing  a  heterogeneous  model  is  to  convert  the  unit  cells  of  the 
lattice  core  to  homogeneous  cells  as  shown  in  Figure  24.  Reactor  cores  are  constructed 
of  several  material  compositions  including  fuel  rods,  cladding,  and  coolant.  Using  the 
general  assumption  that  the  net  neutron  current  flow  across  cell  boundaries  equals  zero, 
one  spatially  averages  the  multigroup  cross  sections  of  the  materials  to  obtain  a  group 
cross  section  for  the  unit  cell.  This  is  usually  done  for  the  fast  and  thermal  group  effects. 
Equation  (54)  defines  the  cell  averaged  group  constant. 
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Hexagonal  Circular  Homogenized 

Assembly  Assembly  Assembly 


Figure  24  Fuel-Cell  Homogenization 


\  dE\  'L{r,E)(l>{r,E)dh 
^  dE  ^  <l>(r,E)d^r 

E,  V„„ 


(54) 


After  homogenizing  the  unit  cells,  control  rods  can  be  added  to  the  homogenized 
core.  One  can  use  the  cell  group  constants  along  with  the  control  rod  cross  sections  in  a 
multigroup  two  dimensional  diffusion  calculation.  These  revised  flux  values  can  then  be 
used  to  calculate  the  final  group  constants  for  the  homogenized  fuel  assembly.  The  final 
step  is  to  calculate  the  flux  and  power  levels  in  the  homogenized  core.  This  can  be 
accomplished  by  dividing  the  core  into  equal  lattice  structures  of  squares  or  other 
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geometric  shapes  that  take  advantage  of  symmetry  to  reduce  the  computation 
requirements. 

The  addition  of  a  thermal-hydraulic  module  would  significantly  enhance  the  model’s 
capability.  PWR  power  distributions  are  coupled  to  temperature.  PWRs  use  water  as  a 
coolant  that  typically  enters  the  bottom  of  the  core  and  leaves  near  the  top.  The  coolant 
decreases  in  density  as  it  absorbs  heat  moving  up  through  the  core.  The  power  density  of 
the  current  model  predicts  a  symmetric  power  peaking  profile  that  is  not  truly  the  case. 
Because  of  the  change  in  density,  the  axial  power  peak  is  actually  slightly  toward  the 
bottom.  The  average  fuel  and  moderator  temperatures  can  then  be  used  to  adjust  the 
macroscopic  cross  sections  for  use  in  the  power  distribution  model. 

Finally,  a  control  adjustment  and  depletion  module  can  be  added.  The  control 
adjustment  module  would  calculate  the  adjustments  necessary  for  control  rod  insertion  or 
withdrawal  to  maintain  criticality  and  the  depletion  module  would  account  for  fuel  bum 
up  impacts  on  the  reactor  core. 
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Appendix  A.  Derivation  of  Three  Dimensional  Source  Integration 


Let 

r  \^^(l)^{r,zydrdGdz  (55) 

for  a  given  homogeneous  material.  Integrating  (p  yields 

/"  =  zydrdz.  (56) 

Let 


/"  (r)  =  {r,z)dz. 


(57) 


Using  the  trapezoid  rule  to  numerically  solve  the  integration  with  respect  to  z,  where  the 
trapezoid  rule  is 


f*b  \r( 


n-1 

s 

/=1 


(58) 


and 


X,  =iAx,  i  =  l,2,3,...n 

A  b 
Ax  =  — 

n 


yields 


r{r)-^2n 


Az 


«>-(r,z„)+r(r,0)+2  ^r{r,z,) 


;=i 


JJ 


(59) 


where 


^o)  “  ®  <^"(^»0)  =  0  due  to  boundary  conditions 

Zj  =  jAz 
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Integrating  with  respect  to  r,  yields 


where 


r 


n-1  ^ 

/"  (0)0  +  /"  (R)R  +  2^ 

,=i  J 


r.  =  iAr 


n 


but  from  equation  (57) 

Zo 

r(R)  =  27ijf{R,z)dz  =  0. 
0 


Simplifying  and  combining  equations  (59)  and  (60)  results  in 

n-1  m-1 

r  =  27rAz(Arf^i^f{r^,Zj) 

i=l  ;=1 


(60) 


(61) 


(62) 
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Appendix  B.  Derivation  of  Right  Circular  Cylinder  Reactor  Core  Solution 


The  basic  equation  is 


D 


(63) 


where  the  boundary  conditions  are 


<Z)(r,±^)  =  0 


Converting  this  into  right  circular  cylinder  coordinates  results  in 


r  dr 


d^) 

_|_ 

rav^ 

dr) 

[Be) 

+  B^^  =  0 


(64) 


where 


Separating  the  variables  and  letting 


<zi(r,z)  =  /?(r)Z(z) 


yields 


(rR'(r)Z(z))+R(r)Z'(z)  +  B^R(r)Z(z)  =  0.  (65) 

r  or 

Collecting  the  terms  provides 


1  (rR')'  Z'(z)  2 

rR(r)  Z(z) 


Now  setting  the  equation  equal  to  a  constant  provides 


(66) 
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(67) 


1  (^^1 

r  R{r) 


+  B^ 


Z{z) 


=  /l‘. 


Solving 


Z'(z)+/l"Z(z)  =  0 


(68) 


where 


yields 


where 


Z 


Z„{z)  =  Cos 


^  nnz^ 

~w 


(69) 


v^y 


n  =  1,3,. 


Now  solving  the  second  equation 


(70) 

r  R{r) 

where 

r[r)=o 

-(rR'^  +  ju^R  =  0  (71) 

r 

and 


The  solution  is  in  the  form  of  the  zeroth  order  Bessel  functions. 
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R{r)  =  AJXtir)  +  CY,{nr) 


(72) 


However  as  r  — >  0,  (//r)  o® ,  therefore  C  must  equal  zero.  At  the  boundary 

condition  r  =  R, 

r(^R^  =  0  =  AJ^(^juR)  (73) 

only  if  A  does  not  equal  zero  and  if  juR  =  ,  where  is  the  zeros  of  the  Jq  . 

Therefore  the  eigenfunctions  and  eigenvalues  are 


for  n  =  0,1,2,... 


(74) 


where  A  is  a  normalization  factor. 
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Appendix  C.  Two  Dimension.  One  Energy  Group  Visual  BASIC  Code 


Thesis  code  by  MAJ  Will  Harman 

This  program  calculates  the  radial  flux/power  profile  for  a  typical  right  circular 
’cylinder  in  two  dimensions.  It  uses  the  standard  diffusion  equation  in  a  one  energy 
group  (one-speed) 

’homogeneous  unreflected  reactor  core.  The  equation  is  solved  by  using  the  finite  central 
’difference  technique.  The  scheme  uses  a  power  iterative  technique  to  solve  for 
’flux  based  upon  an  inital  guess  of  k  effective  and  the  flux.  It  finds  the  eigenfunction 
Tor  the  maximum  eigenvalue  providing  the  fundamental  mode  shape  for  flux. 

Const  mnErrDivByZero  =  11,  mnErrOverFlow  =  6 
Const  mnErrBadCall  =  5 
Private  Sub  Form_Load() 

Toad  initial  values  from  Duderstadt  page  210-211 
Dim  kGuess,  FluxGuess,  SigmaA,  NueSigmaF  As  Double 
Dim  DiffusionCoefficient,  Radius,  h  As  Double 
Textl  =  "" 

Text2  =  0.157  ’NueSigmaF  1/cm 

Texts  =  9.21  DiffusionCoefficient  cm 

Text4  =  120  Radius  cm 

Texts  =  0.5  ’mesh  spacing 

Text6  =  0.9  R  guess 

Texts  =  ""  ’Critical  Rx? 

Text9  =  ""  ’Keffective  will  be  calculated 
TextlO  = " " 

Textl  1  =  "" 

Textl4  = "" 

End  Sub 

Sub  KineticsO 

Dim  prompt  ’ask  user  for  input 
Dim  Valu(2)  As  Double 
Dim  kGuess  As  Double 
Dim  SigmaA  As  Double 

’Set  source  document  for  Excel  chart  to  name  and  location  by  user 
OLE2.SourceDoc  =  ("Textl") 

OLE2.Visible  =  True 
’Ckeck  for  numerical  entries 
ok  =  0 

For  j  =  0  To  1 

ykk  =  Checkin(MaskEdBox(j)) 
ok  =  ok  +  ykk 
Next) 

If  ok  >  0  Then 
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’neutron/cm^2 

1/cm 

1/cm 

’em 

’em 

’spatial  distance  between  nodes  cm 
’number  of  nodes  along  radius 
’em 
’em 

’number  of  iterations  for  convergence  of  k 


’acceptable  error  in  k 

’acceptable  error  in  S  (nueSigmaF*Flux) 


MsgBox  ("You  must  enter  positive  numbers") 

End  If 

’Assign  variables  to  input  boxes 
kGuess  =  MaskEdBox(O) 

FluxGuess  =  MaskEdBox(2) 

SigmaA  =  MaskEdBox(l) 

NueSigmaF  =  CDbl(Text2.Text) 

DiffusionCoefficient  =  CDbl(Text3.Text) 

Radius  =  CDbl(Text4.Text) 
h  =  CDbl(Text5.Text) 
n  =  Radius  /  h 
DeltaR  =  h 
DeltaZ  =  h 
Maxiterations  =  1000 
and  source 
Epsilonk  =  0.00001 
EpsilonS  =  0.015 
kCriticalityTolerance  =  0.0001 

Dim  Sl(),  S2(),  k(),  FluxQ,  EirorSQ,  FluxData()  ’As  Integer 
ReDim  Sl(n) 

ReDim  S2(n) 

ReDim  Flux(n) 

ReDim  k(MaxIterations  +1) 

ReDim  ErrorS(n  - 1) 

ReDim  FluxData(n,  1) 

build  initial  flux  and  source  guess 

Restart  1 :  ’Used  in  ProgramError 

For  i  =  0  To  n  - 1 

Flux(i)  =  FluxGuess 

Sl(i)  =  NueSigmaF  *  Flux(i) 

Next  i 

Restart2:  ’Used  in  ProgramError 
k(0)  =  kGuess 

On  Error  GoTo  ProgramError 

test  =  1  /  k(0)  ’Check  if  user  input  k  value. 

’!  Outer  iterations 

For  i  =  0  To  Maxiterations 

Elux  is  the  only  term  coming  out  of  Call  statement 

Call  LinearFiniteDifference(Sl,  n,  k(i),  SigmaA,  DiffusionCoefficient,  h,  Radius,  Flux) 
’Convert  flux  back  into  source  for  comparison  with  previous  source 
For  m  =  0  To  n 

S2(m)  =  NueSigmaF  *  Flux(m)  ’convert  flux  to  S(n+1;  neutron/cm^3 
Next  m 

’ICalc  area  under  curve  using  composite  trap  rule 
Suml  =  (h  /  2#)  *  S1(0)  ’0 
Sum2  =  (h  /  2#)  *  S2(0)  ’0 
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For  j  =  1  To  n  -  1 
Sum2  =  Sum2  +  h  *  S2(j) 

Suml  =  Suml  +  h  *  Sl(j) 

Next  j 

Integral  S(n+l)/(l/k(n)*Integral  S(n)) 

Build  array  of  Source  errors  to  use  in  tolerance  test 
On  Error  GoTo  ProgramError 

k(i  +  1)  =  k(i)  *  Sum2  /  Suml  Equation  5-274  Duderstadt 
’Calculate  relative  error  between  old  and  new  source 
’!run  to  n-1  because  S=0  at  BC 
For  1  =  0  To  n  - 1 

ErrorS(l)  =  Abs((S2(l)  -  Sl(l))  /  S2(l))  Equation  5-275  Duderstadt 
Next 

Find  maximum  value  of  ErrorS() 

MaxErrorS  =  ErrorS(0) 

For  1  =  1  To  n  - 1 

If  (Errors (1)  >  MaxErrorS)  Then 

MaxErrorS  =  ErrorS(l) 

End  If 
Next  1 

’Check  for  tolerances 

If  ((Abs((k(i  +  1)  -  k(i))  /  k(i))  <  Epsilonk))  And  (MaxErrorS  <  EpsilonS)  Then  Equation 
5-275  Duderstadt 

End  outer  iterations  check  for  convergence  if  true 
kEffective  =  k(i  -i- 1) 

Exit  For 
End  If 

’reset  S1=S2  for  next  iteration 
For  j  =  0  To  n 

sia)=s2a) 

Next  j 

kEffective  =  k(i  -i- 1) 

Next  i  ’end  outer  iteration 

’Check  if  k=l,  if  so  k=keff=critical  Rx 

Numberlterations  =  I 

MultFactor  =  Format(k(i),  "#.#####") 

If  (kEffective  >  1#  -  kCriticalityTolerance)  Then 
Texts  =  "Yes" 

Text9  =  MultFactor  ’k(i) 

Else 

Texts  =  "No" 

Text9  =  MultFactor  ’  k(i) 

’Use  perturbation  to  assist  the  user  in  changing  SigmaA  to  get  criticality 
DeltaSigmaA  =  Abs(l#/  kGuess  - 1#)  *  (-NueSigmaF)  Equation  5-306  Duderstadt 
TextlO  =  SigmaA  +  DeltaSigmaA  ’adjust  new  SigmaA  for  criticality 
End  If 
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If  (i  <  Maxiterations)  Then 
Textl  1  =  Numberlterations 
Else 

Textl4  =  "Exceeded  maxiterations  before  convergence" 

End  If 

Build  chart  using  Excel  and  OLE  capability 
Dim  ExcelApp  As  Object 
Dim  ExcelChart  As  Object 
Dim  ChartTypeVal  As  Integer 

-4100  is  the  value  for  the  MS  Excel  constant  xl3DColumn.  Visual 
BASIC  does  not  understand  MS  Excel  constants,  so  the  value  must  be 
’used  instead. 

’xlLine=4 

’xlXYScatter  =  -4169 
’xl3DSurface=-4103 
Define  my  chart  typ 
ChartTypeVal  =  -4169 

Set  ExcelApp  =  CreateObject("excel. application") 

ExcelApp.Visible  =  False  Will  not  see  Excel  load,  build,  and  chart 
ExcelApp.Workbooks.Add 

’Populate  the  worksheet  in  Excel  with  the  power  (W/cm^3) 

Power  conversion  per  Ott 
For  rwindex  =  0  To  n 

ExcelApp.Cells(rwIndex  +  2,  l).Value  =  h  *  rwindex 

ExcelApp.Cells(rwIndex  +  2,  2). Value  =  Flux(rwlndex)  *  NueSigmaF  /  (2.43  *  3.1  * 
10  ^  10)  W/cm^3 
Next  rwindex 

’select  rows  and  columns  in  worksheet  to  chart  Starts  at  A1  and  highlights  all  values 

ExcelApp.Range("Al").CurrentRegion.Select 

Set  ExcelChart  =  ExcelApp.Charts.Add() 

ExcelChart. Type  =  ChartTypeVal 
ExcelChart.SeriesCollection(l).Name  =  "=""Power""" 

With  ExcelChart 

.Axes(xlCategory,  xlPrimary).HasTitle  =  True 

.Axes(xlCategory,  xlPrimary).AxisTitle.Characters.Text  =  "Radius  (cm)" 
.Axes(xlValue,  xlPrimary).HasTitle  =  True 

.Axes(xlValue,  xlPrimary).AxisTitle.Characters.Text  =  "Power  (Watts/cm'^3)" 

End  With 

ExcelChart.SaveAs  [Textl]  ’Save  chart/data  per  user  input 
’Using  the  square  brackets  tells  Visual  Basic  that  this  is  an 
MS  Excel  command  not  a  Visual  Basic  command. 

OLE2.CreateLink  (Textl)  Link  to  saved  chart 
OLE2. Update  ’allow  immediate  update  of  excel  chart 
ExcelApp.Quit 
Set  ExcelChart  =  Nothing 


Set  ExcelApp  =  Nothing 

’######################################## 

ProgramError: 

Select  Case  Err.Number 
Case  mnErrOverFlow 

MsgBox  ("You  must  guess  an  initial  flux  to  get  a  non-trivial  solution;  the  code  will  guess 
I.OEIO  neutrons/sec-cm^2") 

FluxGuess  =  10000000000# 

Resume  Restart  1 
Case  mnErrDivByZero 

MsgBox  ("You  must  input  a  value  for  k,  the  code  will  assume  kguess=0.9") 

kGuess  =  0.9 

Resume  Restart2 

End  Select 

End  Sub 


Sub  LinearFiniteDifference(Sl,  m,  k,  SigmaA,  DiffusionCoefficient,  h,  Radius,  Flux) 


’ISolves  y(n+l)"=-l/r*y(n+l)’+SigmaA/D*y(n+l)-l/(D*k(n))*NueSigmaF*y(n) 

’!for  y(n+l).  k(n)  and  y(n)  are  calculated  in  main  program. 

’Algorithm  11.3  Burden  &  Faires 
LowerLimit  =  0  ’linner  radius 
UpperLimit  =  m  *  h  ’louter  radius 
alpha  =  0  ’!I.C.  y(LowerLimit)=alpha 

beta  =  0  ’!I.C.  y(upperLimit)=beta 

n  =  m  - 1 

ReDim  A(n)  lower  diagonal 
ReDim  b(n)  ’diagonal 
ReDim  C(n)  ’upper  diagonal 
ReDim  D(n)  ’A.x=d  The  d  vector 
ReDim  l(n) 

ReDim  u(n) 

ReDim  Z(n) 

ReDim  w(n  +1) 

’!Set  distance  of  first  node 
X  =  LowerLimit  +  h  ’em 
’IBuild  diagonals; 

’!a(l)  =  w(i)+w(i-l)  because  flux(0)=flux(l)at  interior  BC 

A(l)  =  2  +  h  ^  2  *  q(x,  SigmaA,  DiffusionCoefficient)  +  (-1  -  (h  /  2#)  *  p(x))  ISlo  units 
b(l)  =  -l  +  (h/  2#)  *  p(x)  ’no  units 

D(l)  =  -h  ^  2  *  r(x,  h,  n,  k,  DiffusionCoefficient,  SI)  ’neutron/cm‘’^2 
For  i  =  2  To  n  - 1 

X  =  LowerLimit  +  i  *  h 

A(i)  =  2  +  h  ^  2  *  q(x,  SigmaA,  DiffusionCoefficient)  ’no  units 
b(i)  =  -1  +  (h  /  2#)  *  p(x)  ’no  units 
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C(i)  =  -1  -  (h  /  2#)  *  p(x)  ’no  units 

D(i)  =  -h  ^  2  *  r(x,  h,  n,  k,  DiffusionCoefficient,  SI)  ’neutrons/cm^2 
Next 

X  =  UpperLimit  -  h 

A(n)  =  2  +  h  ^  2  *  q(x,  SigmaA,  DiffusionCoefficient)  ’no  units 
C(n)  =  -1  -  (h  /  2#)  *  p(x)  ’no  units 

D(n)  =  -h  ^  2  *  r(x,  h,  n,  k,  DiffusionCoefficient,  SI)  +  (1  -  (h  /  2#)  *  p(x))  *  beta 
’neutron/cm^2 

’Crout  Factorization  for  tridiagonal  linear  systems 

’!Back  substitute  for  solution 

1(1)  =  A(l)  ’no  units 

u(l)  =  b(l)  /  A(l)  ’no  units 

Z(l)  =  D(l)  / 1(1)  ’neutron/cm''2 

For  i  =  2  To  n  - 1 

l(i)  =  A(i)  -  C(i)  *  u(i  - 1)  ’no  units 
u(i)  =  b(i)  /  l(i)  ’no  units 

Z(i)  =  (D(i)  -  C(i)  *  Z(i  - 1))  /  l(i)  ’neutron/cm^2 
Next 

l(n)  =  A(n)  -  C(n)  *  u(n  -  1)  ’no  units 

Z(n)  =  (D(n)  -  C(n)  *  Z(n  - 1))  /  l(n)  ’neutron/cm^2 

’!Set  solution  flux  values 

Flux(n  +  1)  =  beta  ’IFlux  at  outer  boundary;  neutron/cm^2 
Flux(n)  =  Z(n)  ’neutron/cm^2 
For  i  =  n  - 1  To  1  Step  -1 
Flux(i)  =  Z(i)  -  u(i)  *  Flux(i  +  1)  ’neutron/cm^2 
Next 

Flux(O)  =  Flux(l)  ’!Set  BC  dFlux/dr=0;  neutron/cm''2 
’End  Subroutine  LinearFiniteDifference 

End  Sub 

’INOTE:  x=radius  of  core  in  cm 
Function  p(x) 

’Real  (8)::  x 
p  =  -1  /  X  ’1/cm 
End  Function 

Function  q(x,  SigmaA,  DiffusionCoefficient) 
q  =  SigmaA  /  DiffusionCoefficient  ’l/cm^2 
End  Function 

Function  r(x,  h,  n,  k,  DiffusionCoefficient,  SI) 

’!S1  is  an  array  filled  by  node  position.  Must  convert  x  to  nodal  points. 

If  (x  =  h)  Then 
X  =  1 

Elself  (x  =  n)  Then 
X  =  n 
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Else 

X  =  X  /  h 

End  If 

r  =  -l/k*Sl(x)*  1/  DiffusionCoefficient  ’neutron/cmM 

End  Function 

Function  Checkin(Box) 

If  Len(Box)  =  0  Then 
Checkin  =  I 
End  If 

End  Function 

Private  Sub  mnu3D2EnergyGroups_Click() 

Load  TwoEnergyGroup 
TwoEnergyGroup.Show 
Unload  k 
End  Sub 

Private  Sub  mnu3DOneEnergyGroupItem_Click() 

Load  rxfrm 
rxfrm.Show 
Unload  k 
End  Sub 

Private  Sub  mnuExitItem_Click() 

End 

End  Sub 

Private  Sub  mnuHelpItem_Click() 

Load  Help 
Help.Show 
End  Sub 

Private  Sub  mnuPrintIteni_Click() 

k.PrintForm 

End  Sub 

Private  Sub  mnuReactorItem_Click() 

Load  Reactor 
Reactor.Show 
End  Sub 

Private  Sub  mnuRunIteni_Click() 

Call  Kinetics 
End  Sub 
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Private  Sub  mnuStartItem_Click() 
Load  ReactorCoreModel 
ReactorCoreModel  .Show 
Unload  k 
End  Sub 

Private  Sub  Option l_Click() 

MsgBox  ("Try  SigmaA=0.1532  1/cm") 
End  Sub 


Appendix  D»  Three  Dimension,  One  Energy  Group  Visual  Basic  Code 


Thesis  code  by  MAJ  Will  Harman 

This  program  calculates  the  radial  and  axial  flux  profile  (3D)  for  a  typical  right 
’circular  cylinder.  It  uses  the  standard  diffusion  equation  in  a  one  energy  group 
’(one-speed)homogeneous  reactor  core.  The  equation  is  solved  by  using  the  finite  central 
’difference  technique  using  a  blocked  tridiagonal  solver.  The  scheme  uses  a  power 
iterative 

’technique  to  solve  for  flux  based  upon  an  inital  guess  of  k  effective  and  the  flux. 

It  finds  the  eigenfunction  for  the  maximum  eigenvalue  providing  the  fundamental  mode 
shape  for  flux. 

’Common  error  statements 
Const  mnSaveAsFailed  =  1004 
Const  mnTypeMismatch  =  13 

Private  Sub  Form_Load() 

’Load  initial  values 

Dim  kGuess,  FluxGuess,  SigmaA,  NueSigmaF  As  Double 

Dim  DiffusionCoefficient,  Radius,  h  As  Double 

Textl  =  180  ’half  Rx  hieght  cm 

Text2  =  0.157  NueSigmaF  1/cm 

Text3  =  9.21  Diffusion  Coefficient  cm 

Text4  =120  ’Core  radius  cm 

Text6  =  "" 

Text?  =  100000000000#  ’  *  10  ^  10  ’neutrons/cm^2 

Texts  =  "" 

Text9  = "" 

TextlO  = "  " 

Textll  =  "" 

Textl4  = "" 

build  fixed  selection  of  mesh  spacing  in  axial  and  radial  cm 
Combo  l.AddItem  "30" 

ComboLAddltem  "20" 

Combo  l.AddItem  "15" 

ComboLAddltem  "12" 

ComboLAddltem  "10" 

Combo  1 . Additem  "6 " 

ComboLAddltem  "5" 

ComboLAddltem  "4" 

ComboLAddltem  "3" 

ComboLAddltem  "2" 

ComboLAddltem  "1" 

Combo2.  Additem  "Choose  half  reactor  core"  ’Select  List  Case  0 
Combo2.  Additem  "Choose  quarter  Rx  core"  ’Select  List  Case  1 
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ComboS.AddItem  "Reduced  Relative  Error;  Slower  Run  Time" 
ComboS.AddItem  "Average  Relative  Error;  Faster  Run  Time" 

End  Sub 

Sub  KineticsQ 
Dim  Valu(2)  As  Double 
Dim  kGuess  As  Double 
Dim  SigmaA  As  Double 
Dim  FluxGuess  As  Double 

’Set  source  document  for  Excel  chart  to  name  and  location  by  user 
OLEl.SourceDoc  =  ("Text6") 

OLEl.  Visible  =  True 
’Ckeck  for  numerical  entries 
ok  =  0 

For  j  =  0  To  1 

ykk  =  Checkin(MaskEdBox(j)) 
ok  =  ok  +  ykk 
Next  j 

If  ok  >  0  Then 

MsgBox  ("You  must  enter  positive  numbers") 

End  If 

’Check  for  ouput  file  location  and  name 
IfText6  =  ""Then 

MsgBox  ("You  must  input  an  output  file  location  and  name") 

End  If 

’Assign  variables  to  input  boxes 
kGuess  =  MaskEdBox(O) 

FluxGuess  =  MaskEdBox(2)  ’CDbl (Text?  .Text)  ’neutron/cm^2 

SigmaA  =  MaskEdBox(l)  ’CDbl(Textl.Text)  ’1/cm 

NueSigmaF  =  CDbl(Text2.Text)  ’1/cm 

DiffusionCoefficient  =  CDbl(Text3.Text)  ’em 

Radius  =  CDbl(Text4.Text)  ’em 

ZHeight  =  CDbl(Textl.Text)  ’em 

On  Error  GoTo  ProgramError 

h  =  CDbl  (Combo  1. Text)  ’spatial  distance  between  nodes;  cm 
n  =  Radius  /  h  ’number  of  nodes  along  radius 
m  =  ZHeight  /  h  ’number  of  nodes  along  z  axis 
DeltaR  =  h  ’em 
DeltaZ  =  h  ’em 

Maxiterations  =  1000  ’number  of  iterations  for  convergence  of  k  and  source 

Provide  the  user  with  a  choice  of  relative  run  times  and  errors 

If  (ComboS.Text  =  "Reduced  Relative  Error;  Slower  Run  Time")  Then 

kTolerance  =  0.0000001 

kCriticalTol  =  0.000001 

Else 

kTolerance  =  0.00001 
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kCriticalTol  =  0.0001 
End  If 

Epsilonk  =  kTolerance  ’acceptable  error  in  k:  Ref  Ott 
EpsilonS  =  0.015  ’acceptable  error  in  S  (nueSigniaF*Flux) 
kCriticalityTolerance  =  kCriticalTol 

Dim  Sl(),  S2(),  k(),  Flux(),  ErrorSQ,  FluxRadialQ,  FluxAxial()  ’As  Integer 
Dim  A(),  1(),  u() 

ReDim  Sl(n  - 1,  m  - 1)  ’S(n) 

ReDim  S2(n  - 1,  m  - 1)  ’S(n+1) 

ReDim  k(MaxIterations  +1) 

ReDim  Flux(n  - 1,  m  - 1) 

ReDim  ErrorS(n  -  1,  m  - 1) 

ReDim  FluxRadial(n,  1) 

ReDim  FluxAxial(m,  1) 

ReDim  A(n  - 1,  m  - 1)  ’stores  main  diagonal  of  matrix 
ReDim  l(n  - 1,  m  - 1)  ’stores  lower  diag  of  tridiag  matrix 
ReDim  u(n  - 1,  m  -  1)  ’stores  upper  diag  of  tridiag  matrix 
’Calculate  initial  source 
For  i  =  0  To  n  -  1 
Forj=0Tom-l  ’(m-l)/2’ 

Flux(i,  j)  =  FluxGuess  ’neutron/cm^2 

Sl(i,  j)  =  NueSigmaF  *  Flux(i,  j)  ’neutron/cm^3; 

Next  j 
Next  i 

’Build  diagonals  of  the  tridiagonal  in  the  blocked  system  by  selecting  reactor 
Select  Case  Combo2.ListIndex  ’either  half  or  quarter  of  Rx.  Boundaries  change. 
Case  0  ’Half  Rx  core 

For  j  =  1  To  m  - 1  ’j=row  position  along  z  axis 
For  i  =  1  To  n  -  1  ’i=column  position  along  radius 
If  i  =  1  Then 

Holdl  =  (-2#  /  DeltaR  ^  2)  -  (2#  /  DeltaZ  ^  2)  -  SigmaA  /  DiffusionCoefficient 
Hold2  =  (1#  /  DeltaR  2)  - 1#  /  (2#  *  DeltaR  2)  ’add  in  boundary  condition 
flux(0,j)=flux(l,j) 

A(i,  j)  =  ((Holdl  +  Hold2)  *  DeltaZ  ^  2)  ’no  units 
Else 

’no  units  for  A(i,j) 

A(i,  j)  =  (((-2#  /  DeltaR  ^  2)  -  2#  /  DeltaZ  ^  2  -  SigmaA  /  DiffusionCoefficient)  * 
DeltaZ  2) 

End  If 
Next  i 
Next) 

Case  1  ’Quarter  Rx  core 
For  j  =  1  To  m  - 1  ’j=row  position  along  z  axis 
For  i  =  1  To  n  - 1  ’i=column  position  along  radius 
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If  (i  =  1  And  j  =  m  - 1)  Then  ’flux(i,j)=flux(i,j-l)=flux(i-l,j) 

Holdl  =  (-2#  /  DeltaR  ^  2)  -  (2#  /  DeltaZ  ^  2)  -  SigmaA  /  DiffusionCoefficient 
Hold2  =  (1#  /  DeltaR  ^  2)  - 1#  /  (2#  *  i  *  DeltaR  ^  2)  ’add  in  boundary  condition 
flux(0,j)=flux(l,j) 

A(i,  j)  =  ((Holdl  +  Hold2  +  1  /  DeltaZ  ^  2)  *  DeltaZ  ^  2)  ’no  units 
Elself  i  =  1  Then  ’flux(i,j)=flux(i-l,j) 

Holdl  =  (-2#  /  DeltaR  ^  2)  -  (2#  /  DeltaZ  ^  2)  -  SigmaA  /  DiffusionCoefficient 
Hold2  =  (1#  /  DeltaR  ^  2)  - 1#  /  (2#  *  i  *  DeltaR  ^  2)  ’add  in  boundary  condition 
flux(0,j)=flux(l,j) 

A(i,  j)  =  ((Holdl  +  Hold2)  *  DeltaZ  ^  2)  ’no  units 
Elself  j  =  m  - 1  Then  ’flux(i,j)=flux(i,j-l) 

A(i,  j)  =  DeltaZ  2  *  ((-2#  /  DeltaR  ^  2)  -  (1#  /  DeltaZ  ^  2)  -  SigmaA  / 
DiffusionCoefficient) 

Else 

A(i,  j)  =  (((-2#  /  DeltaR  ^  2)  -  2#  /  DeltaZ  ^  2  -  SigmaA  /  DiffusionCoefficient)  * 
DeltaZ  ^  2) 

End  If 
Next  i 
Next  j 
End  Select 

’build  lower  diagonal  of  tridiagonal  matrix 
For  j  =  1  To  m  -  1  ’j=row  position  along  z  axis 
For  i  =  1  To  n  -  2  ’i=column  position  along  radius 

l(i,  j)  =  (1#  /  DeltaR  ^  2  - 1#  /  (2#  *  ((i  +  1#)  *  DeltaR  2)))  *  DeltaZ  ^  2  ’no  units 
Next  i 
Next) 

’build  upper  diag.  of  tridiag.  matrix 
For  j  =  1  To  m  - 1 

For  i  =  1  To  n  -  2  ’no  units  for  u(i,j) 

u(i,  j)  =  (1#  /  DeltaR  2  +  1#  /  (2#  *  (i  *  DeltaR  ^  2)))  *  DeltaZ  ^  2 
Next  i 
Next) 

k(0)  =  kGuess 

’[Outer  iterations 

For  i  =  0  To  Maxiterations 

Flux  is  the  only  term  coming  out  of  Call  statement  and  it  is  built 
’so  that  the  first  column  of  the  flux  matrix  equals  the  flux  in  the 
’m-1  row  of  the  Rx  core. 

Call  ThreeDSolver(Sl,  A,  1,  u,  n,  m,  k(i),  SigmaA,  DiffusionCoefficient,  DeltaR,  DeltaZ, 
Flux) 

’build  3D  S2 

Select  Case  Combo2.ListIndex  half  or  quarter  Rx 
Case  0  half  of  Rx 
For  ii  =  0  To  n  - 1 
Forj  =0Tom- 1 
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If  j  =  0  Then 

S2(ii,  0)  =  0  BC 

Elself  (ii  =  0  And  j  <>  0)  Then 

S2(0,  j)  =  NueSigmaF  *  Flux(l,  m  -  j)  ’Flux(0,j)=Flux(l,j) 

Else 

S2(ii,  j)  =  NueSigmaF  *  Flux(ii,  m  -  j) 

End  If 
Next  j 
Next  ii 

Case  I  ’Quarter  of  Rx 

For  ii  =  0  To  n  - 1 

For  j  =  0  To  m  - 1 

If(ii  =  0Andj=0)  Then 

S2(0,  0)  =  NueSigmaF  *  Flux(I,  m  -  1) 

Elself  (ii  =  0  And  j  <>  0)  Then 
S2(0,  j)  =  NueSigmaF  *  Flux(l,  m  -  j) 

Elself  (j  =  0  And  ii  <>  0)  Then 
S2(ii,  0)  =  NueSigmaF  *  Flux(ii,  m  -  1) 

Else 

S2(ii,  j)  =  NueSigmaF  *  Flux(ii,  m  -  j) 

End  If 
Next  j 
Next  ii 
End  Select 

SumI  =  0 
Sum2  =  0 

’Build  3D  integration  of  S(n+I)  and  S(n)  using  composite  trap,  rule 
For  j  =  1  To  m  - 1 
For  b  =  1  To  n  - 1 

Suml  =  Suml  +  b  *  SI(b,  j) 

Sum2  =  Sum2  +  b  *  S2(b,  j) 

Next  b 
Next  j 

’Integral  S(n+l)/(l/k(n)*Integral  S(n))  to  find  next  k  value 
’Build  array  of  Source  errors  to  use  in  tolerance  test 
k(i  +  1)  =  k(i)  *  Sum2  /  Suml  ’Equation  5-275  Duderstadt 
’Calc  3D  relative  error  between  old  and  new  source 
For  ii  =  1  To  n  - 1 
For  j  =  1  To  m  - 1 

ErrorS(ii,  j)  =  Abs((S2(ii,  j)  -  Sl(ii,  j))  /  Abs(S2(ii,  j))) 

Next  j 
Next  ii 

Find  maximum  value  of  ErrorS() 

MaxErrorS  =ErrorS(l,  1) 
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For  ii  =  1  To  n  -  1 
For  j  =  1  To  m  - 1 
If  (ErrorS(ii,  j)  >  MaxErrorS)  Then 
MaxErrorS  =  ErrorS(ii,  j) 

End  If 
Next  j 
Next  ii 

’Check  for  tolerances 

If  ((Abs((k(i  +  1)  -  k(i))  /  k(i))  <  Epsilonk))  And  (MaxErrorS  <  EpsilonS)  Then 
End  outer  iterations  check  for  convergence 
kEffective  =  k(i  +  1) 

Exit  For 
End  If 

’Reassign  S2  to  SI  for  the  next  iteration 
For  ii  =  1  To  m  - 1 
For  j  =  1  To  n  - 1 
Sl(j,ii)  =  S2(j,ii) 

Next  j 
Next  ii 

kEffective  =  k(i  +  1) 

Next  i  ’end  outer  iteration 

’Check  if  k=l,  if  so  k=keff=critical  Rx 

Numberlterations  =  I 

MultFactor  =  Format(k(i),  "#.#####") 

If  (kEffective  >  1#  -  kCriticalityTolerance)  Then 

Texts  =  "Yes" 

Text9  =  MultFactor  ’k(i) 

Else 

Texts  =  "No" 

Text9  =  MultFactor  ’k(i) 

’Use  perturbation  to  assist  the  user  in  changing  SigmaA  to  get  criticality 
DeltaSigmaA  =  Abs(l#  /  kGuess  - 1#)  *  (-NueSigmaF) 

TextlO  =  SigmaA  +  DeltaSigmaA 
End  If 

Textl  1  =  Numberlterations 

’  Textl4  =  "Exceeded  maxiterations  before  convergence" 

Build  Excel  chart  and  spreadsheet 
Dim  ExcelApp  As  Object 
Dim  ExcelChart  As  Object 
Dim  ChartTypeVal  As  Integer 

’-4100  is  the  value  for  the  MS  Excel  constant  xlSDColumn.  Visual 
Basic  does  not  understand  MS  Excel  constants,  so  the  value  must  be 
’used  instead. 

’xlLine=4 
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’xlXYScatter  =  -4169 
’xl3DSurface=-4103 
ChartTypeVal  =  -4103 

Set  ExcelApp  =  CreateObject("excel. application") 

ExcelApp.Visible  =  False  Hide  the  Excel  appliction  from  the  user 
ExcelApp.Workbooks.Add 

Populate  the  Excel  spreadsheet  with  core  power  values  and  locations 
Select  Case  Combo2.ListIndex  Tialf  or  quarter  Rx 
Case  0  Tialf  Rx 
For  rwindex  =  0  To  n 

ExcelApp.Cells(rwIndex  +  2,  l).Value  =  h  *  rwindex 
For  colindex  =  0  To  m 

ExcelApp.Cells(l,  colindex  +  2).Value  =  h  *  colindex 
If  rwindex  =  n  Then 

Excel App.CellsfrwIndex  -i-  2,  colindex  -i-  2). Value  =  0 
Elself  colindex  =  m  Then 

ExcelApp.CellsfrwIndex  +  2,  colindex  -l-  2). Value  =  0 

Elself  colindex  =  0  Then 

Excel App.Cells(rwIndex  +  2,  2). Value  =  0 

Else 

ExcelApp.Cells(rwIndex  -i-  2,  colindex  -i-  2).Value  =  S2(rwlndex,  colindex)  /  (2.43  * 
3.1  *  10  ^  10)  HueSigmaF 
End  If 

Next  colindex 
Next  rwindex 
Case  1  ’Quarter  Rx 
For  rwindex  =  0  To  n 

ExcelApp.Cells(rwIndex  +  2,  l).Value  =  h  *  rwindex 
For  colindex  =  0  To  m 

ExcelApp.Cells(l,  colindex  +  2). Value  =  h  *  colindex 
If  rwindex  =  n  Then 

ExcelApp.Cells(rwIndex  -f-  2,  colindex  +  2). Value  =  0 
Elself  colindex  =  m  Then 

ExcelApp.Cells(rwIndex  +  2,  colindex  +  2).Value  =  0 
Elself  colindex  =  0  Then 

ExcelApp.Cells(rwIndex  +  2,  colindex  +  2).Value  =  S2(rwlndex,  colindex)  /  (2.43 
*3.1*10^10)  NueSigmaF 
Else 

ExcelApp.Cells(rwIndex  -i-  2,  colindex  -i-  2).Value  =  S2(rwlndex,  colindex)  /  (2.43  * 
3.1  *  10  ^  10)  NueSigmaF 
End  If 

Next  colindex 
Next  rwindex 
End  Select 

’select  rows  and  columns  in  worksheet  to  chart 
ExcelApp.Range("Al").CurrentRegion.Select 
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Set  ExcelChart  =  ExcelApp.Charts.Add() 

’Add  legend  information 
ExcelChart.Type  =  ChartTypeVal 
With  ExcelChart 
•HasTitle  =  True 

.ChartTitle.Characters.Text  =  "Power  Plot  of  Reactor  Core" 
.Axes(xlCategory).HasTitle  =  True 

.Axes(xlCategory).AxisTitle.Characters.Text  =  "Height  (cm)" 

.Axes(xlSeries)  .HasTitle  =  True 

.Axes(xlSeries).AxisTitle.Characters.Text  =  "Radius  (cm)" 

.Axes(xlValue).HasTitle  =  True 

.Axes(xlValue).AxisTitle.Characters.Text  =  "Power  (Watts/cm'^'S)" 

End  With 

With  ExcelChart.Axes(xlCategory) 

■HasMajorGridlines  =  False 
•HasMinorGridlines  =  False 
End  With 

With  ExcelChart.Axes(xlSeries) 

•HasMajorGridlines  =  False 
.HasMinorGridlines  =  False 
End  With 

With  ExcelChart.Axes(xlValue) 

.HasMajorGridlines  =  True 
.HasMinorGridlines  =  False 
End  With 

ExcelChart.WallsAndGridlines2D  =  False 
ExcelChart.HasLegend  =  False 
On  Error  GoTo  ProgramError 
ExcelChart.SaveAs  [Text6] 

’Link  Excel  chart  to  saved  name  and  location  for  update 
OLEl.CreateLink  (Text6) 

OLE  1. Update  ’allow  immediate  update  of  excel  chart 
ExcelApp.Quit 

Set  ExcelChart  =  Nothing 
Set  ExcelApp  =  Nothing 

ProgramError: 

Select  Case  Err.Number 
Case  mnTypeMismatch 

MsgBox  ("You  must  select  a  mesh  spacing;  the  code  will  assuem  20  cm") 
h  =  20 

Resume  Next 
Case  mnSaveAsFailed 

MsgBox  ("You  must  provide  a  name  and  file  storage  location  to  update  the  plot  and  store 
data") 
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Resume  Next 
End  Select 
End  Sub 

Function  Checkin(Box) 

If  Len(Box)  =  0  Then 
Checkin  =  1 
End  If 

End  Function 

^  •Im  *14  •!>  «il4  *1^  «kl4  ^  «kl4  vi>  ^1>  •iL*  Vkl*  aX*  •X*  •X*  >1^  vl>  «X*  «X*  "X*  ^  «1>  •!>  ‘X*  ■X*  ^lf**  ^X^  "X*  ^2lf  *X*  ■X^  'X*  ^  S!!>  ^  ^ 

^  ^  ^  rjs  ^|S  ^  ^  ^  ^  ^  «7^  4^  4^  4^  4^*  4|^  4^  4|^  4|^  4^  4jS«  •-J'*  4^  4jV  4jH  4]S  «rj> 

Sub  ThreeDSolver(Sl,  A,  1,  u,  n,  m,  kEffective,  SigmaA,  DiffusionCoefficient,  DeltaR, 
DeltaZ,  X2) 

This  3D  block  tridiagonal  solver  uses  an  iterative  technique  similar  to  Gauss-Seidel. 

It  sets  up  the  banded  tridiagonal  system  based  on  the  discretized  right  circular  cyl 
’diffusion  equation  and  dismantles  that  into  a  blocked  tridiagonal  system. 

Each  block  is  then  solved  with  a  standard  tridiagonal  solver  using  initial  guesses 
’for  the  solution.  The  system  iterates  until  convergence  of  the  solution  (flux  at  the 
’inner  mesh  points) 

ReDim  b(n  - 1,  m  -  1)  ’stores  B  of  A.x=B 

ReDim  X2(n  - 1,  m  - 1)  ’stores  solution  at  mesh  points 

ReDim  XI (n  - 1,  m  - 1)  ’Stores  the  previous  solution  at  mesh  points 

ReDim  e(n  - 1,  m  - 1)  ’stores  rhs  of  block  tridiag  system 

ReDim  AA(n  - 1)  Stores  the  main  diag  (jth  column  of  A  array)  of  block  tridiag 

ReDim  LL(n  - 1)  ’Stores  the  lower  diag  (jth  column  of  1  array)  of  block  tridiag 

ReDim  UU(n  - 1)  ’Stores  the  upper  diag  (jth  column  of  u  array)  of  block  tridiag 

ReDim  EE(n  -  1)  Stores  the  rhs  diag  of  block  tridiag 

ReDim  XX(n)  ’stores  jth  column  solution  vector  from  tridiagonal  solver 

ReDim  ErrorX(n  - 1,  m  -  1)  ’holds  the  max  error  in  convergence  of  solution  in  tridiag 

build  B  vector  of  Ax=B.  B  contains  the  iterative  guess  for  flux. 

’Uses  Sl(i,m-j)  to  convert  the  (i,j)  values  into  m-j  rows  for  the  b  matrix. 

For  j  =  1  To  m  - 1 
For  i  =  1  To  n  -  1 

XI (i,  j)  =  0  ’fill  convergence  test  array  XI  with  0 
X2(i,  j)  =  0  ’fill  with  0 

b(i,  j)  =  DeltaZ  2  *  (-Sl(i,  m  -  j)  /  (kEffective  *  DiffusionCoefficient))  ’neutron/cm'^2 
Next  i 
Next  j 

’solve  block  tridiagonal  system.  See  Solution  Methods  of  my  thesis. 

For  w  =  1  To  10000 
For  i  =  1  To  n  -  1 

e(i,  1)  =  b(i,  1)  -  X2(i,  2)  ’neutron/cm^2 
EE(i)  =  e(i,  1)  ’neutron/cm'^2  typical 

LL(i)  =  l(i,  1)  EL  Lower  diagonal  of  tri diagonal ;no  units  typical 
AA(i)  =  A(i,  1)  ’AA  main  diagonal  of  tridiagoanl;  no  units  typical 
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UU(i)  =  u(i,  1)  TJU  upper  diagonal  of  tridiagonal;  no  units  typical 
Next  i 

Call  Tridiag(n,  LL,  AA,  UU,  EE,  XX)  ’returns  flux  at  row  m-1  in  core 
For  i  =  1  To  n  - 1 

X2(i,  1)  =  XX(i)  ’neutrons/cm^2  typical 
Next! 

For  j  =  2  To  (m  -  2) 

For  i  =  1  To  n  - 1 

e(i,j)  =  b(i,j)-X2(i,j-l)-X2(i,j  +  l) 

EE(i)  =  e(i,j) 

LL(i)  =  l(i,j) 

AA(i)  =  A(i,j) 

UU(i)  =  u(i,j) 

Next  i 

Call  Tridiag(n,  LL,  AA,  UU,  EE,  XX) 

For  i  =  1  To  n  - 1 
X2(i,j)  =  XX(i) 

Next  i 
Nextj 

For  i  =  1  To  n  -  1 

e(i,  m - 1)  =  b(i,  m-1)  - X2(i,  m -  2) 

EE(i)  =  e(i,  m-1) 

LL(i)  =  l(i,  m  - 1) 

AA(i)  =  A(i,  m-1) 

UU(i)  =  u(i,  m  - 1) 

Nexti 

Call  Tiidiag(n,  LL,  AA,  UU,  EE,  XX) 

For  i  =  1  To  n  -  1 
X2(i,  m  - 1)  =  XX(i) 

Nexti 

For  V  =  1  To  n  - 1 
For  j  =  1  To  m  - 1 

EirorXfv,  j)  =  Abs((X2(v,  j)  -  Xl(v,  j))  /  Abs(X2(v,  j))) 

Nextj 
Next  V 

Eind  maximum  value  of  ErrorS() 

MaxErrorX  =  ErrorX(l,  1) 

For  i  =  1  To  n  -  1 

For  j  =  1  To  m  - 1 

Xl(i,j)  =  X2(i,j) ’update  XI 

If  (ErrorX(i,  j)  >  MaxErrorX)  Then 

MaxErrorX  =  ErrorX(i,  j) 

End  If 
Nextj 
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Next  i 

If  (ComboS.Text  =  "Reduced  Relative  Error;  Slower  Run  Time")  Then 

Tolerance  =  0.00000001 

Else 

Tolerance  =  0.001 
End  If 

If  (MaxErrorX  <  Tolerance)  Then 
Exit  For 
End  If 

’Set  boundary  conditions 
Next  w 

If  (w  >  10000)  Then 

MsgBox  ("Exceeded  block  tridiagonal  maxiterations  before  convergence") 
End  If 
End  Sub 


Sub  Tridiag(n,  LL,  AA,  UU,  EE,  XX)  LL=lower  diag,  AA=Diag,  UU=Upper 
diag,EE=A.x 

’Crout  Factorization  for  tridiagonal  linear  systems 
’[Back  substitute  for  solution 
m  =  n-  1 

ReDim  Lower(m) 

ReDim  Upper(m) 

ReDim  Z(m) 

Lower(l)  =  AA(1)  ’no  units 
Upper(l)  =  UU(1)  /  AA(1)  ’no  units 
Z(l)  =  EE(1)  /  Lower(l)  ’neutron/cm^2 
For  i  =  2  To  m  -  1 

Lower(i)  =  AA(i)  -  LL(i  - 1)  *  Upper(i  -  1)  ’no  units 
Upper(i)  =  UU(i)  /  Lower(i)  ’no  units 
Z(i)  =  (EE(i)  -  LL(i  - 1)  *  Z(i  - 1))  /  Lower(i)  ’neutron/cm^2 
Next 

Lower(m)  =  AA(m)  -  LL(m  - 1)  *  Upper(m  - 1)  ’no  units 
Z(m)  =  (EE(m)  -  LL(m  - 1)  *  Z(m  - 1))  /  Lower(m)  ’no  units 
’!Set  solution  flux  values 

XX(m  +  1)  =  0  ’IFlux  at  outer  boundary;  neurton/cm^2 
XX(m)  =  Z(m)  ’neutron/cm^2 
For  i  =  m  - 1  To  1  Step  -1 

XX(i)  =  Z(i)  -  Upper(i)  *  XX(i  +  1)  ’neutron/cm^2 
Next 

XX(0)  =  XX(1)  ’!Set  BC  dFlux/dr=0;  ’neutron/cm^2 
End  Sub  Tridiagonal 
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Private  Sub  mnu2DOneEnergyGroupItem_Click() 

Load  k 

k.Show 

Unload  rxfrm 

End  Sub 

Private  Sub  mnu3DTwoEnergyGroupItem_Click() 

Load  TwoEnergyGroup 

TwoEnergyGroup.Show 

Unload  rxfrm 

End  Sub 

Private  Sub  mnuExitItem_Click() 

End 

End  Sub 

Private  Sub  mnuHelpItem_Click() 

Load  Help 
Help.Show 
End  Sub 

Private  Sub  mnuPrintItem_Click() 

rxfrm.PrintForm 

End  Sub 

Private  Sub  mnuReactorItem_Click() 

Load  Reactor 
Reactor.Show 
End  Sub 

Private  Sub  mnuRunItem_Click() 

Call  Kinetics 
End  Sub 

Private  Sub  mnuStartFormItem_Click() 

Load  ReactorCoreModel 
ReactorCoreModel.Show 
Unload  rxfrm 
End  Sub 

Private  Sub  Option  l_Click() 

MsgBox  ("Try  SigmaA=0.1532  1/cm") 

End  Sub 

Private  Sub  Text5_Change() 
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End  Sub 
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Appendix  E.  Three  Dimension,  Two  Energy  Groups  Visual  Basic  Code 


Thesis  code  by  MAJ  Will  Harman 

This  program  calculates  the  radial  and  axial  Flux/Power  profile  (3D)  for  a  typical  right 
’circular  cylinder  homogeneous  reactor  core  using  two  energy  groups.  It  uses  the 
standard 

’diffusion  equation.  The  equation  is  solved  by  using  the  finite  central 
’difference  technique  using  blocked  tridiagonal  solvwer. 

The  scheme  uses  a  power  iterative  technique  to  solve  for  Fluxl  and  Flux2  based  upon  an 
’inital  guess  of  k  effective  and  Fluxl  and  Flux2.  It  finds  the  eigenfunction  for  the 
maximum 

’eigenvalue  providing  the  fundamental  mode  shape  for  Fluxl  and  Flux2. 


Private  Sub  Form_Load() 

’Load  initial  values 

Dim  kGuess,  FluxlGuess,  SigmaAl,  NueSigmaFl  As  Double 

Dim  DiffusionCoefficientl,  Radius,  h  As  Double 

Textl  =  180  ’half  Rx  hieght  cm 

Text2  =  0.008476  ’NueSigmaFl  1/cm 

Texts  =  1.2627  DiffusionCoefficientl  cm 

Text4  =120  ’Core  radius  cm 

Texts  =  0.18514  ’NueSigmaF2  1/cm 

Text6  =  "" 

Text7  =  0.3543  DiffusionCoefficient2  cm 
Text8  = "" 

Text9  = "" 

TextlO  = "  " 

Textl  1  =  "" 

Text  12  =  0.121  ’SigmaR2  1/cm 
TextlS  =  0.02619  ’SigmaRl 
Textl4  =  0.01207  SigmaAl  1/cm 
TextlS  =  0.121  ’SigmaA2  1/cm 
Textl6  =  0.01412  ’SigmaScatterl2  1/cm 
Textl7  =  0  ’SigmaScatter22  1/cm 

build  fixed  selection  of  mesh  spacing  in  axial  and  radial  cm 
ComboLAddltem  "30" 

Combo  l.AddItem  "20" 

ComboLAddltem  "15" 

ComboLAddltem  "12" 

ComboLAddltem  "10" 

ComboLAddltem  "6" 

ComboLAddltem  "5" 

ComboLAddltem  "4" 

ComboLAddltem  "3" 

ComboLAddltem  "2" 
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Combo  1 . Additem  "  1 " 

Combo2.AddItem  "Choose  half  reactor  core"  ’Select  List  Case  0 
Combo2. Additem  "Choose  quarter  Rx  core"  ’Select  List  Case  1 
ComboS .Additem  "Reduced  Relative  Error;  Slower  Run  Time" 
Combo3. Additem  "Average  Relative  Error;  Faster  Run  Time" 
End  Sub 


Sub  KineticsO 
Dim  Valu(2)  As  Double 
Dim  kGuess  As  Double 
Dim  SigmaAl  As  Double 
Dim  Flux  1  Guess  As  Double 
Dim  Flux2Guess  As  Double 

’Set  source  document  for  Excel  chart  to  name  and  location  by  user 
OLEl.SourceDoc  =  ("Text6") 

OLEl.  Visible  =  True 
OLE2.SourceDoc  =  ("Text6") 

OLE2.Visible  =  True 
’Ckeck  for  numerical  entries 
ok  =  0 

Forj  =0To  1 

ykk  =  Checkin(MaskEdBox(j)) 
ok  =  ok  +  ykk 
Next  j 

If  ok  >  0  Then 

MsgBox  ("You  must  enter  positive  numbers") 

End  If 

’Check  for  ouput  file  location  and  name 
IfText6  =  ""Then 

MsgBox  ("You  must  input  an  output  file  location  and  name") 

End  If 

’Assign  variables  to  input  boxes 
kGuess  =  MaskEdBox(O) 

FluxlGuess  =  MaskEdBox(l)  ’CDbl(Text7.Text)  ’neutron/cm^2 
Flux2Guess  =  MaskEdBox(2) 

SigmaAl  =  CDbl(Textl4.Text)  ’1/cm 

SigmaA2  =  CDbl(Textl5.Text)  ’1/cm 

NueSigmaFl  =  CDbl(Text2.Text)  ’1/cm 

NueSigmaF2  =  CDbl(Text5.Text)  ’1/cm 

DiffusionCoefficientl  =  CDbl(Text3.Text)  ’em 

DiffusionCoefficient2  =  CDbl(Text7.Text)  ’em 

SigmaRl  =  CDbl(Textl3.Text)  ’1/cm  (SigmaTotal-SigmaScatter) 

SigmaR2  =  CDbl(Textl2.Text)  ’1/cm  (SigmaTotal-SigmaScatter) 

SigmaScatterl2  =  CDbl(Textl6.Text)  ’1/cm 

Radius  =  CDbl(Text4.Text)  ’em 
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ZHeight  =  CDbl(Textl.Text)  ’em 

h  =  CDbl(Combol.Text)  ’spatial  distance  between  nodes;  cm 
n  =  Radius  /  h  ’number  of  nodes  along  radius 
m  =  ZHeight  /  h  ’number  of  nodes  along  z  axis 
DeltaR  =  h  ’em 
DeltaZ  =  h  ’em 

Maxiterations  =  1000  ’number  of  iterations  for  convergence  of  k  and  source 

Provide  the  user  a  choice  of  run  times  and  relative  error 

If  (Combo3.Text  =  "Reduced  Relative  Error;  Slower  Run  Time")  Then 

kTolerance  =  0.0000001 

kCriticalTol  =  0.000001 

Else 

kTolerance  =  0.00001 
kCriticalTol  =  0.0001 
End  If 

Epsilonk  =  kTolerance  ’acceptable  error  in  k;  Ref  Ott 
EpsilonS  =  0.015  ’acceptable  error  in  S  (NueSigmaFl*Fluxl) 
kCriticalityTolerance  =  kCriticalTol 

Dim  Sl(),  S2(),  k(),  FluxlQ,  ErrorSQ  ’,  FluxlRadialQ,  FluxlAxialQ  ’As  Integer 
Dim  Flux2(),  Flux2Radial(),  Flux2Axial()  ’As  Integer 
DimAlQ,  A2(),  1(),  u() 

ReDim  Sl(n  - 1,  m  - 1)  ’S(n) 

ReDim  S2(n  - 1,  m  - 1)  ’S(n+1) 

ReDim  k(MaxIterations  +  1) 

ReDim  Fluxl(n  - 1,  m  - 1) 

ReDim  Flux2(n  - 1,  m  - 1) 

ReDim  ErrorS(n  - 1,  m  - 1) 

ReDim  Al(n-l,m-l)  ’stores  main  diagonal  of  matrix 

ReDim  A2(n  - 1,  m  - 1)  ’stores  main  diagonal  of  matrix 

ReDim  l(n  - 1,  m  -  1)  ’stores  lower  diag  of  tridiag  matrix 

ReDim  u(n  - 1,  m  -  1)  ’stores  upper  diag  of  tridiag  matrix 

’Calculate  initial  source  two  D 

For  i  =  0  To  n  - 1 

For  j  =  0  To  m  - 1  ’(m  - 1)  /  2  ’ 

Fluxl(i,  j)  =  FluxlGuess  ’neutron/cm^2 
Flux2(i,  j)  =  Flux2Guess  ’neutron/cm^2 

Sl(i,  j)  =  NueSigmaF2  *  Flux2(i,  j)  +  NueSigmaFl  *  Fluxl(i,  j)  ’neutron/cm^3; 

Cos(3.141592654  *  j  *  DeltaZ  /  Height)  *  Bessel 

Next) 

Next  i 

’Build  diagonals  of  the  tri diagonal  in  the  blocked  system  for  Fluxl 

Select  Case  Combo2.ListIndex  ’either  half  or  quarter  of  Rx.  Boundaries  change. 

Case  0  ’Half  Rx  core 

For  j  =  1  To  m  - 1  ’j=row  position  along  z  axis 
For  i  =  1  To  n  - 1  ’i=column  position  along  radius 
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Ifi  =  lThen 

Holdl  =  (-2#  /  DeltaR  2)  -  (2#  /  DeltaZ  ^  2)  -  SigmaRl  /  DiffusionCoefficientl 
Hold2  =  (1#  /  DeltaR  ^  2)  - 1#  /  (2#  *  i  *  DeltaR  ^  2)  ’add  in  boundary  condition 
Fluxl(0,j)=Fluxl(l,j) 

Al(i,  j)  =  ((Holdl  +  Hold2)  *  DeltaZ  ^  2)  ’no  units 
Else 

’no  units  for  A(i,j) 

Al(i,  j)  =  (((-2#  /  DeltaR  ^  2)  -  2#  /  DeltaZ  2  -  SigmaRl  /  DiffusionCoefficientl)  * 
DeltaZ  ^  2) 

End  If 
Next  i 
Next  j 

Case  1  ’Quarter  Rx  core 
For  j  =  1  To  m  - 1  ’j=row  position  along  z  axis 
For  i  =  1  To  n  -  1  ’i=column  position  along  radius 
If  (i  =  1  And  j  =  m  - 1)  Then  Eluxl(i,j)=Fluxl(i,j-l)=Fluxl(i-l,j) 

Holdl  =  (-2#  /  DeltaR  2)  -  (2#  /  DeltaZ  2)  -  SigmaRl  /  DiffusionCoefficientl 
Hold2  =  (1#  /  DeltaR  ^  2)  - 1#  /  (2#  *  i  *  DeltaR  ^  2)  ’add  in  boundary  condition 
Fluxl(0,j)=Fluxl(l,j) 

Al(i,  j)  =  ((Holdl  +  Hold2  +  1  /  DeltaZ  ^  2)  *  DeltaZ  ^  2)  ’no  units 
Elself  i  =  1  Then  Eluxl(i,j)=Fluxl(i-l,j) 

Holdl  =  (-2#  /  DeltaR  ^  2)  -  (2#  /  DeltaZ  2)  -  SigmaRl  /  DiffusionCoefficientl 
Hold2  =  (1#  /  DeltaR  ^  2)  - 1#  /  (2#  *  i  *  DeltaR ''  2)  ’add  in  boundary  condition 
Fluxl(0,j)=Fluxl(l,j) 

Al(i,  j)  =  ((Holdl  +  Hold2)  *  DeltaZ  ^  2)  ’no  units 
Elself  j  =  m  -  1  Then  Eluxl(i,j)=Fluxl(i,j-l) 

Al(i,  j)  =  DeltaZ  2  *  ((-2#  /  DeltaR  ^  2)  -  (1#  /  DeltaZ  ^  2)  -  SigmaRl  / 
DiffusionCoefficient  1 ) 

Else 

Al(i,  j)  =  (((-2#  /  DeltaR  ^  2)  -  2#  /  DeltaZ  2  -  SigmaRl  /  DiffusionCoefficientl)  * 
DeltaZ  2) 

End  If 
Next  i 
Next  j 
End  Select 

’Build  diagonal  of  the  tridiagonal  for  Flux2 

Select  Case  Combo2.ListIndex  ’either  half  or  quarter  of  Rx.  Boundaries  change. 

Case  0  ’Half  Rx  core 

For  j  =  1  To  m  - 1  ’j=row  position  along  z  axis 
For  i  =  1  To  n  - 1  I=column  position  along  radius 
Ifi  =  1  Then 

Holdl  =  (-2#  /  DeltaR  2)  -  (2#  /  DeltaZ  2)  -  SigmaR2  /  DiffusionCoefficient2 
Hold2  =  (1#  /  DeltaR  ^  2)  - 1#  /  (2#  *  i  *  DeltaR  ^  2)  ’add  in  boundary  condition 
Fluxl(0,j)=Fluxl(l,j) 

A2(i,  j)  =  ((Holdl  +  Hold2)  *  DeltaZ  ^  2)  ’no  units 
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Else 

’no  units  for  A(i,j) 

A2(i,  j)  =  (((-2#  /  DeltaR  ^  2)  -  2#  /  DeltaZ  ^  2  -  SigmaR2  /  DiffusionCoefficient2)  * 
DeltaZ  ^  2) 

End  If 
Next  i 
Next  j 

Case  1  ’Quarter  Rx  core 
For  j  =  1  To  m  - 1  ’j=row  position  along  z  axis 
For  i  =  1  To  n  -  1  ’i=column  position  along  radius 
If  (i  =  1  And  j  =  m  - 1)  Then  FluxI(i,j)=FIuxl(i,j-I)=FluxI(i-l,j) 

Holdl  =  (-2#  /  DeltaR  ^  2)  -  (2#  /  DeltaZ  ^  2)  -  SigmaR2  /  DiffusionCoefficient2 
Hold2  =  (1#  /  DeltaR  ^  2)  - 1#  /  (2#  *  i  *  DeltaR  ^  2)  ’add  in  boundary  condition 
Fluxl(OJ)=Fluxl(l,j) 

A2(i,  j)  =  ((Holdl  +  Hold2  +  1  /  DeltaZ  ^  2)  *  DeltaZ  ^  2)  ’no  units 
Elself  i  =  1  Then  Eluxl(i,j)=Fluxl(i-l,j) 

Holdl  =  (-2#  /  DeltaR  ^  2)  -  (2#  /  DeltaZ  ^  2)  -  SigmaR2  /  DiffusionCoefficient2 
Hold2  =  (1#  /  DeltaR  2)  - 1#  /  (2#  *  i  *  DeltaR  ^  2)  ’add  in  boundary  condition 
Fluxl(0,j)=Fluxl(l,j) 

A2(i,  j)  =  ((Holdl  +  Hold2)  *  DeltaZ  ^  2)  ’no  units 
Elself  j  =  m  - 1  Then  Fluxl(i,j)=Fluxl(i,j-l) 

A2(i,  j)  =  DeltaZ  ^  2  *  ((-2#  /  DeltaR  ^  2)  -  (1#  /  DeltaZ  2)  -  SigmaR2  / 
DiffusionCoefficient2) 

Else 

A2(i,  j)  =  (((-2#  /  DeltaR  2)  -  2#  /  DeltaZ  ^  2  -  SigniaR2  /  DiffusionCoefficient2)  * 
DeltaZ  ^  2) 

End  If 
Next  i 
Next  j 
End  Select 

’build  lower  diagonal  of  tridiagonal  matrix 
For  j  =  1  To  m  - 1  ’j=row  position  along  z  axis 
For  i  =  1  To  n  -  2  ’i=colutnn  position  along  radius 

l(i,  j)  =  (1#  /  DeltaR  ^  2  - 1#  /  (2#  *  ((i  +  1#)  *  DeltaR  ^  2)))  *  DeltaZ  ^  2  ’no  units 
Nexti 
Next  j 

’build  upper  diag.  of  tridiag.  matrix 
For  j  =  1  To  m  -  1 

For  i  =  1  To  n  -  2  ’no  units  for  u(i,j) 

u(i,  j)  =  (1#  /  DeltaR  ^  2  +  1#  /  (2#  *  (i  *  DeltaR  ^  2)))  *  DeltaZ  ^  2 
Next  i 
Next  j 

k(0)  =  kGuess 
’!  Outer  iterations 
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For  i  =  0  To  Maxiterations 

Fluxl  is  the  only  term  coming  out  of  Call  statement  and  it  is  built 
’so  that  the  first  column  of  the  Fluxl  matrix  equals  the  Fluxl  in  the 
’m-1  row  of  the  Rx  core. 

SolvingFlux  =  1  ’Selection  of  B  vector  in  A.x=B 
’Solve  for  fast  flux  values 

Call  ThreeDSolver(SolvingFlux,  SI,  Al,  1,  u,  n,  m,  k(i),  DeltaR,  DeltaZ,  Fluxl,  Fluxl) 
SolvingFlux  =  2  ’Selection  of  B  vector  in  A.x=B 
’Use  fast  flux  values  and  solve  for  thermal  values 

Call  ThreeDSolver(SolvingFlux,  SI,  A2, 1,  u,  n,  m,  k(i),  DeltaR,  DeltaZ,  Fluxl,  Flux2) 
’build  3D  S2 

Select  Case  Combo2.ListIndex  half  or  quarter  Rx 

Case  0  half  of  Rx 

For  ii  =  0  To  n  - 1 

For  j  =  0  To  m  - 1 

Ifj  =  0Then 

S2(ii,  0)  =  0  BC 

Elself  ii  =  0  Then 

S2(0,  j)  =  NueSigmaFl  *  Fluxl(l,  m  - 1)  +  NueSigmaF2  *  Flux2(l,  m  -  j) 

Flux  l/2(0,j)=Flux  1/2(1, j) 

Else 

S2(ii,  j)  =  NueSigmaFl  *  Fluxl(ii,  m  -  j)  +  NueSigmaF2  *  Flux2(ii,  m  -  j) 

End  If 
Next] 

Next  ii 

Case  1  ’Quarter  of  Rx 
For  ii  =  0  To  n  - 1 
For  j  =  0  To  m  - 1 
If(ii  =  0Andj=0)Then 

S2(0,  0)  =  NueSigmaFl  *  Fluxl(l,  m - 1)  +  NueSigmaF2  * Flux2(l,  m-1) 

Elself  ii  =  0  Then 

S2(0,  j)  =  NueSigmaFl  *  Fluxl(l,  m  -  j)  +  NueSigmaF2  *  Hux2(l,  m  -  j) 

Elself  (j  =  0  And  ii  <>  0)  Then 

S2(ii,  0)  =  NueSigmaFl  *  Fluxl(ii,  m  -  1)  +  NueSigmaF2  *  Flux2(ii,  m-1) 

Else 

S2(ii,  j)  =  NueSigmaFl  *  Fluxl(ii,  m  -  j)  +  NueSigmaF2  *  Flux2(ii,  m  -  j) 

End  If 
Next] 

Next  ii 
End  Select 

Suml  =  0 
Sum2  =  0 

Build  3D  integration  of  S(n+1)  and  S(n)  using  composite  trap,  rule 
For  j  =  1  To  m  - 1 
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For  b  =  1  To  n  - 1 

Suml  =  Suml  +  b  *  Sl(b,  j) 

Sum2  =  Sum2  +  b*S2(b,j) 

Next  b 
Next  j 

Integral  S(n+l)/(l/k(n)*Integral  S(n))  to  find  next  k  value 
Build  array  of  Source  errors  to  use  in  tolerance  test 
k(i  +  1)  =  k(i)  *  Sum2  /  Suml  Equation  5-275  Duderstadt 
’Calc  3D  relative  error  between  old  and  new  source 
For  ii  =  1  To  n  - 1 
For  j  =  1  To  m  - 1 

ErrorS(ii,  j)  =  Abs((S2(ii,  j)  -  Sl(ii,  j))  /  Abs(S2(ii,  j))) 

Next  j 
Next  ii 

Find  maximum  value  of  ErrorS() 

MaxErrorS  =  ErrorS(l,  1) 

For  ii  =  1  To  n  -  1 

For  j  =  1  To  m  - 1 

If  (ErrorS(ii,  j)  >  MaxErrorS)  Then 

MaxErrorS  =  ErrorS(ii,  j) 

End  If 
Next  j 
Next  ii 

’Check  for  tolerances 

If  ((Abs((k(i  +  1)  -  k(i))  /  k(i))  <  Epsilonk))  And  (MaxErrorS  <  EpsilonS)  Then 
End  outer  iterations  check  for  convergence 
kEffective  =  k(i  +  I) 

Exit  For 
End  If 

’Reassign  S2  to  SI  for  the  next  iteration 
For  ii  =  1  To  m  - 1 
For  j  =  1  To  n  - 1 
Sl(j,ii)  =  S2(j,ii) 

Next  j 
Next  ii 

kEffective  =  k(i  +  I) 

Next  i  ’end  outer  iteration 

’Check  if  k=I,  if  so  k=keff=critical  Rx 

Numberlterations  =  I 

MultFactor  =  Format(k(i),  "#.#####") 

If  (kEffective  >  I#  -  kCriticalityTolerance)  Then 
Texts  =  "Yes" 

Text9  =  MultFactor  ’k(i) 

Else 

Texts  =  "No" 

Text9  =  MultFactor  ’  k(i) 


S9 


MsgBox  ("The  system  is  not  critical.  Please  try  changing  the  core  composition  density" 
&_ 

"  or  core  geometry") 

End  If 

Text  11  =  Numberlterations 
If  (i>  1000)  Then 

MsgBox  ("Exceeded  maxiterations  before  convergence") 

End  If 

’######################################### 

Build  Excel  chart  and  spreadsheets 
Dim  ExcelApp  As  Object 
Dim  ExcelChartl  As  Object 
Dim  ExcelChart2  As  Object 
Dim  ChartTypeVal  As  Integer 

-4100  is  the  value  for  the  MS  Excel  constant  xlSDColumn.  Visual 
Basic  does  not  understand  MS  Excel  constants,  so  the  value  must  be 
’used  instead. 

’xlLine=4 

’xlXYScatter  =  -4169 
’xl3DSurface=-4103 
ChartTypeVal  =  -4169  ’-4103 
Set  ExcelApp  =  CreateObject("excel.application") 

ExcelApp.Visible  =  False 
ExcelApp.Workbooks.Add 

’Allow  the  user  to  choose  which  node  to  plot  the  data  on  both  radial  and  axial 
Dim  prompt  1,  prompt2 

promptl  =  "The  number  of  interior  radial  mesh  spaces  ="  &  n  _ 

&  ".  Please  choose  the  mesh  point  between  0  and  "  &  n  &  "  to  plot  the  axial  power." 
plotAxial  =  InputBox$(promptl) 

prompt2  =  "The  number  of  interior  axial  mesh  spaces  ="  &  m  _ 

&  ".  Please  choose  the  interior  mesh  point  between  0  and  "  &  m  &  "  to  plot  the  radial 
power." 

plotRadial  =  InputBox$(prompt2) 

Select  Case  Combo2.ListIndex  Iialf  or  quarter  Rx 

Case  0  ’half  Rx 

For  rwindex  =  0  To  n 

Bill  Excel  sheet  1  with  Radial  power  data 

ExcelApp.Sheets("Sheetl").Cells(rwIndex  +  2,  1). Value  =  h  *  rwindex 
ExcelApp.Sheets("Sheetl").Cells(rwIndex  +  2,  2). Value  =  Flux  1  (rwindex,  m  - 
plotRadial)  *  NueSigmaFl  /  (2.43  *  3.1  *  10  ^  10)  ’W/cm''3 

ExcelApp.Sheets("Sheetl").Cells(rwIndex  +  2,  3).Value  =  Flux2(rwlndex,  m  - 
plotRadial)  *  NueSigmaF2  /  (2.43  *  3.1  *  10  ^  10)  ’W/cm^3 

ExcelApp.Sheets("Sheetl").Cells(rwIndex  +  2,  4). Value  =  (Flux  1  (rwindex,  m  - 
plotRadial)  *  NueSigmaFl  +  Flux2(rwlndex,  m  -  plotRadial)  *  NueSigmaF2)  /  (2.43  * 
3.1  *  10  ^  10)  ’W/cm^3 
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Next  rwindex 

For  colindex  =  0  To  m 

’fill  Excel  sheet2  with  axial  power  data 

ExcelApp.Sheets("Sheet2").Cells(m  -  colindex  +  2,  l).Value  =  h  *  colindex 
ExcelApp.Sheets("Sheet2").Cells(colIndex  +  2,  2).Value  =  Fluxl(plotAxial, 
colindex)  *  NueSigmaFl  /  (2.43  *  3.1  *  10  ^  10)  ’W/cm^3 

ExcelApp.Sheets("Sheet2").Cells(colIndex  +  2,  3).Value  =  Flux2(plotAxial, 
colindex)  *  NueSigmaF2  /  (2.43  *  3.1  *  10  10)  ’W/cm^3 

ExcelApp.Sheets("Sheet2").Cells(colIndex  +  2, 4).Value  =  (Flux  1  (plot  Axial, 
colindex)  *  NueSigmaFl  +  Flux2(plotAxial,  colindex)  *  NueSigmaF2)  /  (2.43  *  3.1  *  10 
^  10)  ’W/cm^3 
Next  colindex 
Case  1  ’Quarter  Rx 
For  rwindex  =  0  To  n 

ExcelApp.Sheets("Sheetl").Cells(rwIndex  +  2, 1). Value  =  h  *  rwindex 
ExcelApp.Sheets("Sheetl").Cells(rwIndex  +  2,  2).Value  =  Flux  1  (rwindex,  m  - 
plotRadial)  *  NueSigmaFl  /  (2.43  *  3.1  *  10  ^  10)  ’W/cm^3 

ExcelApp.Sheets("Sheetl").Cells(rwIndex  +  2,  3). Value  =  Flux2(rwlndex,  m  - 
plotRadial)  *  NueSigmaF2  /  (2.43  *  3.1  *  10  ^  10)  ’W/cm^3 

ExcelApp.Sheets("Sheetl’').Cells(rwIndex  +  2, 4).Value  =  (Flux  1  (rwindex,  m  - 
plotRadial)  *  NueSigmaFl  +  Flux2(rwlndex,  m  -  plotRadial)  *  NueSigmaF2)  /  (2.43  * 
3.1  *  10  ^  10)  ’W/cm^3 
Next  rwindex 
For  colindex  =  0  To  m 

ExcelApp.Sheets("Sheet2").Cells(m  -  colindex  +  2,  l).Value  =  h  *  colindex 
ExcelApp.Sheets("Sheet2").Cells(colIndex  +  2,  2).Value  =  Flux  1  (plot Axial, 
colindex)  *  NueSigmaFl  /  (2.43  *  3.1  *  10  ^  10)  ’W/cm'^3 

ExcelApp.Sheets("Sheet2").Cells(eolIndex  +  2,  3).Value  =  Flux2(plotAxial, 
colindex)  *  NueSigmaF2  /  (2.43  *  3.1  *  10  10)  ’W/cm^3 

ExcelApp.Sheets("Sheet2").Cells(colIndex  +  2,  4).Value  =  (Flux  1  (plot Axial, 
colindex)  *  NueSigmaFl  +  Flux2(plot Axial,  colindex)  *  NueSigmaF2)  /  (2.43  *  3.1  *  10 
^  10)  ’W/cm^3 
Next  colindex 
End  Select 

’select  rows  and  columns  in  worksheet  to  chart 

ExcelApp.Sheets("Sheetl").Range("Al").CurrentRegion.Select 

Set  ExcelChartl  =  Excel  App.Charts.Add() 

ExcelApp.Sheets("Sheet2").Range("Al").CurrentRegion.Select 

ExcelChartl. Type  =  ChartTypeVal 

ExcelChartl. SeriesCollection(l).Name  =  "=""Thermal’""' 

ExcelChartl. SeriesCollection(2).Name  =  "=""Fast""" 
ExcelChartl.SeriesCollection(3).Name  =  "=""Totar""' 

’  ExcelChart.Location  Where:=xlLocationAsObject,  Name:="Sheetl" 

With  ExcelChartl 
.HasTitle  =  True 

.ChartTitle.Characters.Text  =  "Power  in  Reactor  Core" 
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.Axes(xlCategory,  xlPrimary).HasTitle  =  True 

.Axes(xlCategory,  xlPrimary).AxisTitle.Characters.Text  =  "Radius  (cm)" 
.Axes(xlValue,  xlPrimary).HasTitle  =  True 

.Axes(xl  Value,  xlPrimary).AxisTitle.Characters.Text  =  "Power  (Watts/cm^3)" 

End  With 

ExcelChartl.HasLegend  =  True 
Build  Chart  2  for  Core  Height  profile 

ExcelApp.Sheets("Sheet2").Select  ’.Range("Ar’).CurrentRegion. Select 
ExcelApp.Sheets("Sheet2").Range("Al").CurrentRegion. Select 
Set  ExcelChart2  =  Excel  App.Charts.Add() 

ExcelChart2.Type  =  ChartTypeVal 
ExcelChart2.SeriesCollection(l).Name  =  "=""Thermal""" 
ExcelChart2.SeriesCollection(2).Name  =  "=""Fast""" 
ExcelChart2.SeriesCollection(3).Name  =  "=""Total""" 

With  ExcelChart2 
.HasTitle  =  True 

•ChartTitle.Characters.Text  =  "Power  in  Reactor  Core" 

.Axes(xlCategory,  xlPrimary).HasTitle  =  True 

.Axes(xlCategory,  xlPrimary).AxisTitle.Characters.Text  =  "Core  Height(cm)" 
.Axes(xlValue,  xlPrimary).HasTitle  =  True 

.Axes(xl Value,  xlPrimary).AxisTitle.Characters.Text  =  "Power  (Watts/cm^3)" 

End  With 

ExcelChart2.HasLegend  =  True 

’save  chart,  activate  chart,  OLE  link  chart,  OLE  update  chart  for  chart  1  &  2 
ExcelChartl.SaveAs  [Text6] 

ExcelApp.Sheets("Chartl").Select  ’activate  chart 
OLEl.CreateLink  (Text6) 

OLEl. Update  ’allow  immediate  update  of  excel  chart 
ExcelApp.Sheets("Chart2").Select  ’Activate  chart 
OLE2.CreateLink  (Text6) 

OLE2.Update 
ExcelApp.Quit 
Set  ExcelChart  =  Nothing 
Set  ExcelApp  =  Nothing 

End  Sub 

Function  Checkin(Box) 

If  Len(Box)  =  0  Then 
Checkin  =  1 
End  If 

End  Function 

Sub  ThreeDSolver(SolvingFlux,  SI,  A,  1,  u,  n,  m,  kEffective,  DeltaR,  DeltaZ,  Fluxl,  X2) 
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This  3D  block  tridiagonal  solver  uses  an  iterative  technique  similar  to  Gauss-Seidel. 
It  sets  up  the  banded  tridiagonal  system  based  on  the  discretized  right  circular  cyl 
’diffusion  equation  and  dismantles  that  into  a  blocked  tridiagonal  system. 

’Each  block  is  then  solved  with  a  standard  tridiagonal  solver  using  initial  guesses 
’for  the  solution.  The  system  iterates  until  convergence  of  the  solution  (Fluxl  at  the 
Inner  mesh  points) 

SigmaScatterl2  =  CDbl(Textl6.Text)  ’1/cm 
DiffusionCoefficientl  =  CDbl(Text3.Text)  ’em 
DiffusionCoefficient2  =  CDbl(Text7.Text)  ’em 
ReDim  b(n  - 1,  m  - 1)  ’stores  B  of  A.x=B 
ReDim  X2(n,  m)  ’stores  solution  at  mesh  points 
ReDim  Xl(n,  m)  ’Stores  the  previous  solution  at  mesh  points 
ReDim  e(n  -  1,  m  - 1)  ’stores  rhs  of  block  tridiag  system 
ReDim  AA(n  - 1)  ’Stores  the  main  diag  (jth  column  of  A  array)  of  block  tridiag 
ReDim  LL(n  - 1)  ’Stores  the  lower  diag  (jth  column  of  1  array)  of  block  tridiag 
ReDim  UU(n  -  1)  ’Stores  the  upper  diag  (jth  column  of  u  array)  of  block  tridiag 
ReDim  EE(n  - 1)  ’Stores  the  rhs  diag  of  block  tridiag 
ReDim  XX(n)  ’stores  jth  column  solution  vector  from  tridiagonal  solver 
ReDim  ErrorX(n  -  1,  m  -  1)  ’holds  the  max  error  in  convergence  of  solution  in  tridiag 
’Fill  solutions  with  0 
For  j  =  0  To  m 
For  i  =  0  To  n 

Xl(i,  j)  =  0  ’fill  convergence  test  array  XI  with  0 
X2(i,  j)  =  0 ’fill  with  0 
Next  i 
Next  j 

’build  B  vector  of  Ax=B.  B  contains  the  iterative  guess  for  Fluxl. 

If  SolvingFlux  =  1  Then 
For  j  =  1  To  m  - 1 
For  i  =  1  To  n  -  1 

b(i,  j)  =  DeltaZ  ^  2  *  (-Sl(i,  m  -  j)  /  (kEffective  *  DiffusionCoefficientl)) 
’neutron/cm^2 
Next  i 
Next  j 

Else  ’SolvingFlux  =  2 
For  j  =  1  To  m  - 1 
For  i  =  1  To  n  - 1 

b(i,  j)  =  DeltaZ  2  *  (-SigmaScatterl2  *  Fluxl(i,  j)  /  DiffusionCoefficient2) 
’neutron/cm^2 
Next  i 
Next  j 
End  If 

’solve  block  tridiagonal  system.  See  Solution  Methods  in  my  thesis. 

For  w  =  1  To  100000 
For  i  =  1  To  n  -  1 

e(i,  1)  =  b(i,  1)  -  X2(i,  2)  ’neutron/cm^2 
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EE(i)  =  e(i,  1)  ’neutron/cm'^2  typical 

LL(i)  =  l(i,  1)  ’LL  Lower  diagonal  of  tridiagonal;no  units  typical 
AA(i)  =  A(i,  1)  ’AA  main  diagonal  of  tridiagoanl;  no  units  typical 
UU(i)  =  u(i,  1)  ’UU  upper  diagonal  of  tridiagonal;  no  units  typical 
Next! 

Call  Tridiag(n,  LL,  AA,  UU,  EE,  XX)  ’returns  Fluxl  at  row  m-1  in  core 
For  i  =  1  To  n  - 1 

X2(i,  1)  =  XX(i)  ’neutrons/cm^2  typical 
Next  i 

For  j  =  2  To  (m  -  2) 

For  i  =  1  To  n  - 1 

e(i,j)  =  b(i,j)-X2(i,j-l)-X2(i,j  +  l) 

EE(i)  =  e(i,j) 

LL(i)  =  l(i,j) 

AA(i)  =  A(i,j) 

UU(i)  =  u(i,j) 

Next  i 

Call  Tridiag(n,  LL,  AA,  UU,  EE,  XX) 

For  i  =  1  To  n  -  1 
X2(i,j)  =  XX(i) 

Next  i 
Next) 

For  i  =  1  To  n  - 1 

e(i, m - 1)  =  b(i,  m-1) - X2(i,  m -  2) 

EE(i)  =  e(i,  m  - 1) 

LL(i)  =  l(i,  m-1) 

AA(i)  =  A(i,  m-1) 

UU(i)  =  u(i,  m  - 1) 

Next! 

Call  Tridiag(n,  LL,  AA,  UU,  EE,  XX) 

For  i  =  1  To  n  - 1 
X2(i,  m  -  1)  =  XX(i) 

Next  i 

For  i  =  0  To  m  -  1 
X2(0,  i)  =  X2(l,i) 

Next  i 

For  i  =  0  To  n 

Select  Case  Combo2.ListIndex  ’either  half  or  quarter  of  Rx.  Boundaries  change. 

Case  0  ’Half  Rx  core 

X2(i,  m)  =  0  ’update  XI 

Case  1  ’Quarter  Rx 

X2(i,  m)  =  X2(i,  m  - 1)  ’update  XI 

End  Select 

X2(i,  0)  =  0 

Next! 
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For  V  =  1  To  n  -  1 
For  j  =  1  To  m  - 1 

ErrorX(v,  j)  =  Abs((X2(v,  j)  -  Xl(v,  j))  /  Abs(X2(v,  j))) 

Next  j 
Next  V 

Find  maximum  value  of  ErrorS() 

MaxErrorX  =  ErrorX(l,  1) 

For  i  =  1  To  n  - 1 

For  j  =  1  To  m  - 1 

Xl(i,j)  =  X2(i,j)  ’update  XI 

If  (ErrorX(i,  j)  >  MaxErrorX)  Then 

MaxErrorX  =  ErrorX(i,  j) 

End  If 
Next  j 
Next  i 

If  (ComboS.Text  =  "Reduced  Relative  Error;  Slower  Run  Time")  Then 

Tolerance  =  0.00000001 

Else 

Tolerance  =  0.001 
End  If 

If  (MaxErrorX  <  Tolerance)  Then 
Exit  For 
End  If 

’Set  boundary  conditions 
Next  w 

If  (w>  10000)  Then 

MsgBox  ("Exceeded  block  tridiagonal  maxiterations  before  convergence") 
End  If 
End  Sub 


Sub  Tridiag(n,  LL,  AA,  UU,  EE,  XX)  ’LL=lower  diag,  AA=Diag,  UU=Upper 
diag,EE=A.x 

’Crout  Factorization  for  tridiagonal  linear  systems 
’IBack  substitute  for  solution 
m  =  n  - 1 
ReDim  Lower(m) 

ReDim  Upper(m) 

ReDim  Z(m) 

Lower(l)  =  AA(1)  ’no  units 
Upper(l)  =  UU(1)  /  AA(1)  ’no  units 
Z(l)  =  EE(1)  /  Lower(l)  ’neutron/cm^2 
For  i  =  2  To  m  - 1 

Lower(i)  =  AA(i)  -  LL(i  - 1)  *  Upper(i  - 1)  ’no  units 
Upper(i)  =  UU(i)  /  Lower(i)  ’no  units 
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Z(i)  =  (EE(i)  -  LL(i  - 1)  *  Z(i  - 1))  /  Lower(i)  ’neutron/cin^2 
Next 

Lower(m)  =  AA(m)  -  LL(m  - 1)  *  Upper(m  - 1)  ’no  units 
Z(m)  =  (EE(m)  -  LL(m  - 1)  *  Z(m  - 1))  /  Lower(m)  ’no  units 
’!Set  solution  Fluxl  values 

XX(m  +  1)  =  0  ’IFluxl  at  outer  boundary;  neurton/cm^2 
XX(m)  =  Z(m)  ’neutron/cm''2 
For  i  =  m  - 1  To  1  Step  -1 
XX(i)  =  Z(i)  -  Upper(i)  *  XX(i  +  1)  ’neutron/cm^2 
Next 

XX(0)  =  XX(1)  ’!Set  BC  dFluxl/dr=0;  ’neutron/cm^2 
End  Sub  Tridiagonal 

Private  Sub  Option  l_Click() 

MsgBox  ("Try  SigmaAl=0.1532  1/cm") 

End  Sub 

Private  Sub  mnu2DOneEnergyGroupItem_Click() 

Load  k 
k.Show 

Unload  TwoEnergyGroup 
End  Sub 

Private  Sub  mnu3DOneEnergyGroupItem_Click() 

Load  rxfrm 
rxfrm.Show 

Unload  TwoEnergyGroup 
End  Sub 

Private  Sub  mnuExitItem_Click() 

End 

End  Sub 

Private  Sub  mnuHelpItem_Click() 

Load  Flelp 
Help.Show 
End  Sub 

Private  Sub  mnuPrintItem_Click() 
TwoEnergyGroup.PrintForm 
End  Sub 

Private  Sub  mnuReactorItem_Click() 

Load  Reactor 
Reactor.Show 
End  Sub 
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Private  Sub  mnuRunItem_Click() 

Call  Kinetics 
End  Sub 

Private  Sub  mnuStari;FormItem_Click() 

Load  ReactorCoreModel 
ReactorCoreModel  .Show 
Unload  TwoEnergyGroup 
End  Sub 

Private  Sub  OLEl_Updated(Code  As  Integer) 
OLEl.  Visible  =  True 
End  Sub 

Private  Sub  OLE2_Updated(Code  As  Integer) 
OLE2.Visible  =  True 
End  Sub 
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Appendix  F.  Sample  Output  Charts  and  Data 

The  figures  below  are  typical  charts  provided  as  output  to  the  user.  In  addition  to  the 
charts,  the  program  saves  the  data  on  an  Excel  worksheet  for  later  use  by  the  program 
user. 
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Figure  25  Two  Dimensional  Output  with  Mesh  Spacing  =  0.5  cm 
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Power  Plot  of  Reactor  Core 
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Figure  26  3D,  One  Energy  Group,  Half  Core  Plot  with  Mesh  Spacing  =  6  cm 

Power  Plot  of  Reactor  Core 
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Figure  28  3D,  Output  with  Mesh  Spacing  =  6  cm,  Two  Energy  Groups 
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Core  Height(cm) 


Figure  29  3D,  Output  with  Mesh  Spacing  =  6  cm,  Two  Energy  Groups 
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Figure  30  3D,  Output  with  Mesh  Spacing  =  6  cm,  Two  Energy  Groups  Half  Rx 


Power  in  Reactor  Core 
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Figure  31  3D,  Output  with  Mesh  Spacing  =  6  cm.  Two  Energy  Groups  Half  Rx 


Appendix  G.  Relative  Error  Plots  of  Test  Cases 


3D,  1  Energy  Group,  Quarter  Rx 
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Figure  32  Radial  Plot,  Mesh  Spacing  =  30  cm 


3D,  1  Energy  Group,  Quarter  Rx 


Numerical  Solution 
Mathematical  Solution 


Figure  33  Axial  Plot,  Mesh  Spacing  =  30  cm 
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Normalized  Value  I  Normalized  Value 


3D,  1  Energy  Group,  Quarter  Rx 
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Figure  34  Radial  Plot,  Mesh  Spacing  =  15  cm 


Figure  35  Axial  Plot,  Mesh  Spacing  =  15  cm 
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3D,  1  Energy  Group,  Quarter  Rx 


Radius  (cm) 


Figure  36  Radial  Plot,  Mesh  Spacing  =  10  cm 


3D,  1  Energy  Group,  Quarter  Rx 


Figure  37  Axial  Plot,  Mesh  Spacing  =  10  cm 
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3D,  1  Energy  Group,  Quarter  Rx 


Radius  (cm) 


Figure  38  Radial  Plot,  Mesh  Spacing  =  5  cm 


Figure  39  Axial  Plot,  Mesh  Spacing  =  5  cm 
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Normalized  Values  Normalized  Value 


Figure  42  Radial  Plot,  Mesh  Spacing  =  3  cm 


3D,  1  Energy  Group,  Quarter  Rx 


Figure  43  Axial  Plot,  Mesh  Spacing  =  3  cm 


3D,  1  Energy  Group,  Quarter  Rx 


Radius  (cm) 


Figure  44  Radial  Plot,  Mesh  Spacing  =  2  cm 


3D,  1  Energy  Group,  Quarter  Rx 


Height  (cm) 


Figure  45  Axial  Plot,  Mesh  Spacing  =  2  cm 
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Normalized  Value 


Figure  46  Radial  Plot,  Mesh  Spacing  =  1  cm 


3D,  1  Energy  Group,  Quarter  Rx 


Figure  47  Axial  Plot,  Mesh  Spacing  =  1  cm 


19 


3D,  Total  2  Energy  Groups,  Half  Rx 


Figure  48  Radial  Plot,  Mesh  Spacing  =  30  cm 


Figure  49  Axial  Plot,  Mesh  Spacing  =  30 
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3D,  Total  2  Energy  Group,  Half  Rx 


— Numerical  Solution 

— Mathematical 
Solution 


Figure  50  Radial  Plot,  Mesh  Spacing  =  20  cm 


3D,  Total  2  Energy  Groups,  Half  Rx 


— Numerical  Solution 

— Mathematical 
Solution 


Figure  51  Axial  Plot,  Mesh  Spacing  =20 
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Figure  52  Radial  Plot,  Mesh  Spacing  =  15  cm 


Figure  53  Axial  Plot,  Mesh  Spacing  =  15  cm 
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Normalized  Value 


Figure  54  Radial  Plot,  Mesh  Spacing  =  10  cm 


3D,  Total  2  Energy  Groups,  Half  Rx 


Height  (cm) 


Figure  55  Axial  Plot,  Mesh  Spacing  =  10  cm 
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Figure  56  Radial  Plot,  Mesh  Spacing  =  5  cm 


3D,  Total  2  Energy  Groups,  Half  Rx 
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Figure  57  Axial  Plot,  Mesh  Spacing  =  5  cm 
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Figure  58  Radial  Plot,  Mesh  Spacing  =  4  cm 


Figure  59  Axial  Plot,  Mesh  Spacing  =  4  cm 


3D,  Total  2  Energy  Groups,  Half  Rx 


Height  (cm) 


Figure  61  Axial  Plot,  Mesh  Spacing  =  3  cm 
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Figure  62  Radial  Plot,  Mesh  Spacing  =  2  cm 


3D,  Total  2  Energy  Groups,  Half  Rx 


Figure  63  Axial  Plot,  Mesh  Spacing  =  2  cm 


17 


3D,  Total  2  Energy  Groups,  Half  Rx 


Radius  (cm) 


Figure  64  Radial  Plot,  Mesh  Spacing  =  1  cm 


3D,  Total  2  Energy  Groups,  Half  Rx 


Figure  65  Axial  Plot,  Mesh  Spacing  =  1  cm 
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Appendix  H.  Linking  Visual  BASIC  and  Excel 


There  are  several  references  available  to  assist  in  creating  embedded  OLE  Excel 
charts  within  Visual  BASIC;  however,  they  do  not  fully  explain  how  to  chart  an  array  of 
data.  This  appendix  provides  the  general  procedure  to  chart  and  embed  an  OLE  Excel 
object  consisting  of  an  array  of  data  as  well  provide  tips  to  creating  more  complicated 
charts  and  links. 

The  key  to  creating  an  embedded  OLE  Excel  chart  is  to  understand  the  OLE 
commands  and  the  Excel  application  commands.  To  set  the  OLE  source  documentation 
to  link  the  OLE  to  the  source  file,  provide  the  name  and  location  of  the  Excel  file  that 
will  contain  the  data  generated  by  the  Visual  BASIC  code. 

OLE.SourceDoc  =  (“Text6”) 

In  this  example,  “Text6”  is  the  TextBox  on  the  Visual  BASIC  form  that  the  program  user 
uses  to  input  the  name  and  location  of  the  Excel  file.  SourceDoc  is  a  procedure  that 
links  the  OLE  to  the  source  document.  To  make  the  OLE  object  visible  on  the  Visual 
BASIC  Form,  use  the  following  code. 

OLE.Visible  =  True 

This  should  be  placed  in  the  code  so  that  the  OLE  object  becomes  visible  only  when 

desired.  Visual  BASIC  is  object  oriented.  Objects  have  built-in  procedures  and  settings 

that  allow  the  programmer  to  control  the  functionality  of  the  object.  To  build  the  Excel 

chart,  dimension  each  object  as  shown  below. 

Dim  ExcelApp  As  Object 
Dim  ExcelChart  As  Object 

This  will  allow  each  of  the  newly  defined  objects  to  have  an  associated  property  or 
method  drop  down  window  displaying  the  available  commands  in  Visual  BASIC. 
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To  define  the  Excel  chart  type,  Visual  BASIC  must  be  given  the  Excel  constants  instead 
of  the  chart  name.  Some  common  Excel  constants  are  given  in  Table  12. 


Table  12  Excel  Constants  for  Charts 


Chart  Type 

Excel  Constant 

xlLine 

4 

xlColumn 

3 

xlXYScatter 

-4169 

xlSDBar 

-4099 

xlSDSurface 

-4103 

To  create  a  three  dimensional  chart,  dimension  the  variable  name  for  the  chart  type 

and  assign  it  an  Excel  constant  value. 

Dim  ChartTypeVal  As  Integer 
ChartTypeVal  =  -4103 

Build  an  Excel  Workbook  and  Worksheet  using  the  following  commands. 

Set  ExcelApp  =  CreateObject("excel. application") 

ExcelApp.Visible  =  False 
ExcelApp.Workbooks.Add 

This  adds  a  workbook  to  Excel  and  keeps  the  Excel  code  running  in  the  background 
without  being  visible  to  the  program  user.  To  see  Excel  run  during  the  Visual  BASIC 
runtime  mode,  change  "false"  to  "true".  This  can  assist  the  programmer  during 
debugging  because  it  allows  the  program  user  to  see  how  the  data  is  being  added  the 
worksheet. 

To  add  data  to  the  Excel  Worksheet,  use  a  For-Next  loop  as  shown. 
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For  rwlndex  =  0  to  n 

ExcelApp.Sheets("Sheet1  ").Cells(rwlndex,collndex).Value  =  your  data 
Next  rwindex 

This  adds  the  data  to  Worksheet  one  in  the  cells  corresponding  to  the  Excel  (row, 
column)  coordinate  system.  Once  the  data  is  added  to  the  Worksheet  it  can  then  be 
charted  and  linked  to  the  Visual  BASIC  OLE.  The  following  sample  code  selects  the 
data  on  the  sheet,  defines  the  chart  type,  adds  data  series  and  axis  labels,  saves  the  chart, 
and  links  the  chart  to  the  OLE. 

ExcelApp.Sheets("Sheet1  ").Range("A1  ").CurrentRegion.Select 
Set  ExcelChart!  =  ExcelApp.Charts.AddQ 
ExcelChartI  .Type  =  ChartTypeVal 
ExcelChart1.SeriesCollection(1).Name  =  "="  "Thermal""" 

ExcelChartI  .SeriesCollection(2). Name  =  "=""Fast . 

ExcelChartI  .SeriesCollection(3).Name  =  "=""Total . 

With  ExcelChartI 
.HasTitle  =  True 

■ChartTitle.Characters.Text  =  "Power  in  Reactor  Core" 

.Axes{xlCategory,  xIPrimary). HasTitle  =  True 

.Axes(xlCategory,  xlPrimary).AxisTitle.Characters.Text  =  "Radius  (cm)" 
.Axes(xlValue,  xIPrimary). HasTitle  =  True 

.Axes(xlValue,  xlPrimary).AxisTitle.Characters.Text  =  "Power  (Watts/cm^3)" 
End  With 

ExcelChartI  .HasLegend  =  True 

’save  chart,  activate  chart,  OLE  link  chart,  OLE  update  chart  for  chart  1 
ExcelChartI  .SaveAs  [Texts] 

ExcelApp.Sheets("Chart1  ").Select  ’activates  the  desired  chart 

OLE1  .CreateLink  (Texts)  ’Creates  the  link  to  the  name  and  location 

’given  in  TextBox  six 

OLE1  .Update  ’allow  immediate  update  of  excel  chart  on 

’the  Visual  BASIC  Form 

ExcelApp.Quit  ’Quits  the  application 

Set  ExcelChart  =  Nothing  ’  Clears  the  previous  settings 

Set  ExcelApp  =  Nothing 

This  sample  code  produced  the  OLE  embedded  object  shown  in  Figure  66  directly  on 
the  Visual  BASIC  Form.  To  open  Excel  and  access  the  chart  and  data,  double  click  the 
OLE  on  the  Visual  BASIC  Form. 
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Power  Plot  of  Reactor  Core 


Radius  (cm> 


Figure  66  Sample  OLE  Embedded  Object 

Excel  contains  Visual  BASIC  with  Applications  that  is  a  limited  version  of  Visual 
BASIC.  Although  the  Visual  BASIC  with  Applications  will  not  run  many  Visual  BASIC 
commands,  it  can  assist  it  developing  Visual  Basic  code  needed  to  create  embedded 
charts.  For  example,  to  develop  sample  Visual  BASIC  codes  use  the  macro  command 
while  in  Excel  to  record  the  steps  in  building  an  Excel  chart  and  then  display  the  code 
using  Visual  BASIC  with  Applications.  This  will  provide  the  general  coding  language  to 
develop  variations  to  Excel  charts.  In  some  cases,  the  sample  code  displayed  in  Visual 
BASIC  with  Applications  can  be  copied  directly  into  Visual  BASIC. 
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